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Walsh functions form an orthonormal basis set consisting of square waves. The dis- 
continuous nature of square waves make the system well suited for representing functions 
with discontinuities. The product of any two Walsh functions is another Walsh function — 
a feature that can radically change an algorithm for solving non-linear partial differential 
equations (PDEs). The solution algorithm of non-linear differential equations using Walsh 
function series is unique in that integrals and derivatives may be computed using simple ma- 
trix multiplication of series representations of functions. Solutions to PDEs are derived as 
functions of wave component amplitude. Three sample problems are presented to illustrate 
the Walsh function series approach to solving unsteady PDEs. These include an advection 
equation, a Burgers equation, and a Riemann problem. The sample problems demonstrate 
the use of the Walsh function solution algorithms, exploiting Fast Walsh Transforms in 
multi-dimensions (0(iVlog(iV))). Details of a Fast Walsh Reciprocal, defined here for the 
first time, enable inversion of a Walsh Symmetric Matrix in 0(iVlog(iV)) operations. Walsh 
functions have been derived using a fractal recursion algorithm and these fractal patterns 
are observed in the progression of pairs of wave number amplitudes in the solutions. These 
patterns are most easily observed in a remapping defined as a fractal fingerprint (FFP). A 
prolongation of existing solutions to the next highest order exploits these patterns. The 
algorithms presented here are considered a work in progress that provide new alternatives 
and new insights into the solution of non-linear PDEs. 


I. Introduction 

The solution of partial differential equations (PDEs) with Walsh functions offers new opportunities to 
simulate many challenging problems in mathematical physics. The approach was developed to better simulate 
hypersonic flows with shocks on unstructured grids. Fundamental derivations of Walsh functions, their 
products, reciprocals, derivatives, and integrals have been documented . 1 Examples of shocked flows in 
nozzles were used to demonstrate associated algorithms. 

While this algorithm holds great promise, there are coding infrastructure challenges that make it difficult 
to further explore and develop numerical simulation tools using Walsh functions. Some unique strengths of 
Walsh functions are the capabilities to systematically account for non-linear interactions of component waves 
and to explicitly propagate boundary conditions across the solution domain. However, the implementation 
of these capabilities presents tedious programming challenges. Multiplication and division of two Walsh 
functions are fundamentally different from currently supported operations in any programming language. 
Keeping track of Jacobians required for the solution of systems of PDEs can become especially tedious in 
this environment. Fortunately, all of this tedious detail can be supported in the background through use of a 
module that supports operator overloading and other commonly used functions and routines in the solution 
of partial differential equations. This support greatly simplifies the learning curve for new users wanting to 
explore the use of Walsh functions in solving PDEs. A FORTRAN module for supporting Walsh function 
simulations has been documented . 2 

The purpose of this paper is to document algorithm advances associated with the solution of unsteady 
PDEs, with emphasis on the Riemann problem. A longer term goal of this work is to evaluate a multi- 
dimensional Riemann solver with small basis set to replace existing flux reconstruction algorithms on tetra- 
hedral grids. In addition, two supporting algorithms not previously documented are provided within: 
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• The Fast Walsh Transform (FWT) in multi-dimensions is key to efficient problem solving using a large 
number ( N ) of terms in the series. An element of any solution algorithm requires the transformation 
from wave number representation to discrete value representation and back. The FWT implements 
this transformation in order N\og(N ) operations (if N = 2 P ). 

• The Fast Walsh Reciprocal (FWR) provides the inverse of a Walsh symmetric matrix in order Nlog(N) 
operations. Any problems involving the evaluation of the reciprocal of a function requires evaluation 
of a Walsh symmetric matrix inverse. 

This paper is organized as follows. Section II provides a review of Walsh function fundamentals. This 
review is based on a more comprehensive article 1 published in February, 2014. Section III presents algorithms 
for filtering high frequency noise and wave component sequencing (equivalent to grid sequencing in a finite- 
difference context). Section IV defines the three unsteady test problems and documents how FORTRAN 
module WALSH_TOOLS 2 is used to obtain their solutions. Both an implicit and explicit formulation are 
described. The implicit formulation uses an exact Jacobian of the Walsh series components. It converges 
very well but becomes intractably large for more than 2 13 degrees of freedom in a sub-domain. The explicit 
formulation utilizes an approximate Jacobian with respect to the full series representation on a sub-domain. 
Critical algorithms affecting solution times are described in the Appendices. The Fast Walsh Transform 
(FWT) in multi-dimensions is described in Sec. VII (Appendix A). The Fast Walsh Reciprocal (FWR) is 
described in Sec. VIII (Appendix B). 


II. Walsh Function Fundamentals 


A. Walsh Function Derivation 


Let g n (x) be a basis function over the domain x a < x < Xb with n segments. The orthonormal basis function 
requirements 3 that 

rxb 

/ g n {x)g m (x)dx = S nm (1) 

J X a 

enable a recursive algorithm to define g n (x) based on the following constraints. Constrain \g n (x)\ equal to 
a constant across the entire domain. Let index p > 0 define a segment size dx p . 


dXp — {Xb Xa)/2 P 


(2) 


An additional constraint on segment lengths within g n (x) requires that at most two segment lengths ( dx p and 
2 dx p ) are allowed to create n segments spanning the domain. These new constraints allow basis functions 
to be grouped according to the smallest segment size. A basis function g n (x) belongs to group p if 


2 P ~ 1 

As an example, consider the evolution of the first 
four groups of g n (x) corresponding to 0 < p < 3 on 
the domain with x a = 0 and Xb = 1. For p = 0 
the only element of the group according to Eq. 3 is 
n = 1. A single segment with length given by Eq. 

2 spans the domain. Therefore, gi(x) = 1 (Fig. 1) 
is the first basis function, belonging to group p = 

0 and satisfying Eq. 1. (While g \(x) = — 1 also 
satisfies Eq. 1, it is convenient to adopt a convention 
where the sign of the function on the first segment 
is positive.) 

For p = 1, the only element of the group accord- 
ing to Eq. 3 is n = 2. Two segments with length equal 


<n<2 p (3) 

g,(x) 
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Figure 1. Basis function in group p= 0. 
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where function gi (x) (Fig. 2) is the second basis function, belonging to group p = 1 and satisfying Eq. 1. 
The value of the function at interior segment boundaries is set to 0, the average of the function on either 
side of the boundary. (The vertical line through interior segment boundaries in Fig. 2 provides strong visual 
impact of the discontinuity and is not meant to imply all values between -1 and 1 are included.) 

The group defined by p = 2 is the first to have 
more than one function, gs(x) and g^x). The al- g 2 ( x ) 


lowed segment sizes in this group are 1/4 and 1/2. 1 

Each function is initialized with the smallest seg- 
ment at the beginning and end of the domain. For 0 



n = 3, this initialization leaves a single segment 
of length 1/2 available to span the interior of the 
domain. For n = 4, the initialization leaves two -1 

segments of length 1/4 to span the interior of the 

) 0.25 0 

5 x 0.75 i 


domain. The corresponding functions are shown in 

Fig. 3. By inspection, it is clear that these first four Figure 2. Basis function in group p=l. 

basis functions all satisfy Eq. 1. 
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Figure 3. Basis functions in group p= 2. 


The group defined by p = 3 includes four functions, g§(x) through g% (a:). The allowed segment sizes in 
this group are 1/8 and 1/4. Each function is again initialized with the smallest segment at the beginning 
and end of the domain. This is the first group where the distribution of segments required to satisfy Eq. 1 
for all combinations of functions in groups 0-3 becomes non-trivial. Still, the number of permutations of 
segment size distribution is small and a solution, evident by inspection, is shown in Fig. 4. 

These results suggest the existence of an algorithm that would allow the definition of basis functions 
for p » 3 and, correspondingly n » 8. Such an algorithm was first identified by Walsh 3 4 and the basis 
functions are used extensively in signal and image processing. 5,6 A new derivation applies a fractal-like 
bisection algorithm focused on the distribution of segment sizes that yields the following definition for basis 
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Figure 4. Basis functions in group p=3. 


functions g n (x) on the domain x a < x < x^'. 


( — l)(m+l)( Xb _ Xa y 1/2 

for 

'Em— 1 y X <C. Xjyi 

0 < m < n 

0 

for 

x = Xm 

0 < m < n 

(• Xb - Xg )^ 1 / 2 

for 

X = X a 


k (~ 1)(" +1 ) (X b - X a )~ 1/2 

for 

X = Xb 



where 

m 

x m = x a + y ^g n {i)dxp 
1=1 


( 4 ) 


( 5 ) 


and i, to, and n are positive integers. In Eq. 5, g n (i) is a vector of length n in which the value of element i 
signifies the length of the i th segment in g n (x). Thus, it defines a dimensionless distribution such that 


/ ■% I 1 if x i Xi — i dxp 

9n{i) = < where 1 < i < n (6) 

I 2 if Xi — Xi- 1 = 2 dx p 

and dx p was defined previously in Eq. 2. 

The dimensionless distribution function g n {i) is defined as a function of group number p. The sum of the 
segment lengths must equal the domain length, consequently it follows that 

n 

y, g n [i) = ( X b - Xa)/dXp = 2 P (7) 

i=l 

Figure 5 presents the sequence of the first 32 g n segment length distribution vectors in a manner that 
highlights the symmetries inherent in the defining algorithm. (Details of the folding algorithm to obtain 
this result are provided in the original reference. 1 ) Note that the two columns on the right side of the 
figure include the index n of the basis function and its group number p. Recall that the minimum segment 
length in a basis function dx p decreases by a factor of 2 for each increment in group number p. Compare the 


4 of 44 


American Institute of Aeronautics and Astronautics 







n p 

1 1 0 

1 1 2 1 

12 1 3 2 

1111 4 2 

1 2 2 2 1 5 3 

121121 6 3 

1112111 7 3 

1111111 1 8 3 

12222222 1 9 4 

1 2 2 2 1 1 2 2 2 1 10 4 

1211222112 1 11 4 

1 2 1 1 2 1 1 2 1 1 2 1 12 4 

1112112112111 13 4 

11121111112111 14 4 

111111121111111 15 4 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 16 4 

12222222222222221 17 5 

122222221122222221 18 5 

1222112222222112221 19 5 

12221122211222112221 20 5 

12112221122211222112 1 21 5 

1211222112112112221121 22 5 

12112112112221121121121 23 5 

121121121121121121121121 24 5 

1112112112112112112112111 25 5 

11121121121111112112112111 26 5 

111211111121121121111112111 27 5 

1112111111211111121111112111 28 5 

11111112111111211111121111111 29 5 

111111121111111111111121111111... 30 5 

1111111111111112111111111111111.. 31 5 

11111111111111111111111111111111. 32 5 


Figure 5. The distribution vector g n for level p ranging from 0 to 5 and 1 < rt < 2 s . 


distribution of segments in Figs. 1 - 4 to the top eight rows of Fig. 5 to confirm understanding of the relation 
between the basis function g n (x) and the corresponding segment distribution vector g n {i). Figure 6 shows 
the segments used to partition the interval x a < x < Xb for the first 64 basis functions g n (x). The symbols 
on each line indicate segment boundaries. The crease lines for the first three bisecting folds are indicated 
along the bottom of the figure. A more careful examination of Fig. 6 reveals that the folding algorithm used 
to define the distribution of segments propagates patterns through ever smaller bisections of the domain and 
larger values of n. For example, take a sheet of paper to cover the right side of Fig. 6 from bisection 1. The 
pattern that remains for even (2m) or odd (2m-l) values of n on the bisected figure is exactly that of the 
whole figure. The only difference between these same two even and odd basis functions on the full domain is 
that the even function had two segments of length dx p on either side of the crease (bisection 1) and the odd 
function had a single segment of length 2dx p spanning the crease. Take another sheet of paper and cover 
the left side of the figure from bisection 2. The pattern that remains for every fourth value of n is exactly 
that of the whole. These patterns have been exploited 1 to derive a new recursion formula for g n (x) that 
lead to proofs for closure under multiplication and explain the relationships of specific pairs of wave number 
amplitudes to be discussed in a subsequent section. 

B. Product of Two Walsh Functions 

A multiplication table of basis functions is developed in Table 1. The table is constructed based on a 
unit domain where Xb — x a = 1- The appropriate scale factor for the product in the most general case 
is ( Xb — Xa)" 1 / 2 derived from the definition of g n (x) in Eq. 4. This scale factor is inserted in the upper 
right corner of the table. The first row of results is simply a recognition that a constant value across the 
entire domain (g 1 ( 2 ;)) multiplying any other function ( g n {x )) returns the rescaled value of g n [x). The lower 
diagonal of the table indicates that the square of g n {x) on every segment is a constant across the domain 
(except at a finite number of points where g n { x) = 0 on segment boundaries), returning a rescaled value of 
9i{x). 
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Figure 6. Segments spanning domain for the first 64 basis functions g n (x). 


Every other element of Table 1 can be defined by recursion. It is not generally required to perform more 
than one bisection and one union to fill out the table. 1 

C. Series Expansions of Functions 

A function /( x) may be represented by an infinite series 

OO 

^2^ n g n (x) ( 8 ) 

n—1 
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Table 1. A multiplication table of basis functions. 



over the interval x a < x < Xb- The coefficients A n are computed 

r x b 


An — 


f(x)g n (x)dx 


(9) 


Because g n (x) is a constant on any given segment and changes sign across segments, the integral in Eq. 9 
may be expressed as 

n f x i + 1 

An = {x b - x 0 ) _1/2 ^(-l) ( * +1) / f(x)dx (10) 

i = 1 

where 


Xi- f-1 — X a -f“ ^ ( [jn (j)dXp ( 11 ) 

i= i 

Xi — Xi - {-i g n {i)dx p (12) 

Note that for domains comprised of N = 2 P segments the FWT (see Appendix A) can be used to compute 
the wave components A n from the integrated average value of / over each segment in order pN operations. 

Figure 7 shows the series representation of the function f(x) = x on the interval — 1 < x < 1. The basis 
function series is truncated at n = 64 (p = 6). Stair-stepping of the segmented representation is evident; 
nevertheless, the error norm L\ = 0 where Lf is given by 




2 P 

(■ Xb - Xa )- 1 ^ 
m = 1 


2 P 

f(xm+ 1 / 2 ) ^ Aigjjxj 77 ,-|- 1 / 2)1 

i = 1 


dx p 


(13) 


The error norm sums the absolute difference between the function and its series representation at the midpoint 
of the smallest segment size for each segment. Each contribution is weighted by the ratio of dx p to the total 
interval length. A zero norm in this case means that the function and its series representation exactly match 
at each midpoint of the smallest segment in the series. Note that the only non-zero coefficients in Fig. 7b 
are Ak where k = 2,4,8, 16, and 64. All linear functions f{x) = C\X + Co on all intervals may be represented 
using coefficients A ^ with k = 2 P and p > 0 with L\ — 0. 
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(a) f(x) compared to Y?n=i A n9n{x) 


Figure 7. Representation of f(x) 



x by series A n 9 n (x). 




(a) fix) compared to Y^=i A n9n(x) (b) Series coefficients 

Figure 8. Representation of f(x) = cos(a;) by series X]n=i A n9 n ( x )- 


Figure 8 shows the series representation of the function /( x) = cos(x) on the interval — 7r < x < n. 
The error norm for this case is L\ = 0.000255735. With successive increases in p (doubling the number of 
segments) the L\ norm drops by a factor of 4 ( L\ = 0.00006392036, L\ = 0.00001597925). 

Figure 9 shows the series representation of the function f(x) = l — H(x—l/3) on the interval — 1 < x < 1. 
The discontinuity in the Heaviside function is placed at x = 1/3 so that it never falls perfectly on a segment 
boundary. The error norm for this case is L\ = 0.005208338. With successive increases in p (doubling 
number of segments) the Li norm drops by a factor of 2 in this case with a discontinuity ( L\ = 0.002604162, 
L\ = 0.001302088). The stair-stepping across the discontinuity is hardly evident on the full scale using p = 8 
(n = 256 segments), as seen in Fig. 10. There are no overshoots or ringing in the series representation. An 
important detail here is to accurately approximate the integral where the discontinuity crosses a segment. 
The variation of the coefficients A n in Figs. 9(b) and 10(b) show a constant magnitude within a family p 
and a rapid drop in magnitude as p increases. 
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(a) f(x) compared to Y?n=i A n9n{x) 



Figure 9. Representation of f(x) = 1 — H(x - 1/3) by series ^® =1 (*)■ 



(a) f(x) compared to 

Figure 10. Representation of f(x) = 1 — H(x 



- 1/3) by series A n9 n (x). 


D. Series Expansion of Function Reciprocal 

If a function f{x) is represented by the truncated series 

2 P 

22 A n g n (x ) (14) 

n—1 

over the interval x a < x < Xb and its reciprocal, r(x) = l/f(x), is represented by the series 

2 P 

22 B mg m (x) (15) 

m—1 

on the same interval then the coefficients Bm can be computed as a function of coefficients A n as follows: 

B-m ra)] *^a) (Id) 
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where 



A-p( 1 , 1 ) 

A-p( 1,2) 

Av{ 1,3) 

Ap( lj2 p) 


A-p( 2,1) 

A-p( 2,2) 

Av{ 2,3) 

A-p(2, 2P) 

[-^■P(n,ra)] 

Av( 3,1) 

Av( 3,2) 

Ap( 3,3) 

A-p( 3,2P) 


. Ap( 2 P,l) 

A'P(2p,2) 

A'P(2p,3) 

• • • A-p(2P ,2P) 


and V(n, m) refers to the multiplication table shown originally in Table 1 with the mapping explicitly shown 
in Table 4 of Sec. VIII (Appendix B). Note that the inverse of this Walsh symmetric matrix can be computed 
in 0(./V log(iV)) operations as shown in Appendix B. 

As an example, the series expansion for r( x) = 1/x is presented in Fig. 11 over the interval —0.1 < x < 0.2 
as derived from the expansion for f(x ) = x using Eq. 16. 



Figure 11. Representation of r(x ) 



(b) p = 10, 1024 segments 

1/x derived from the series expansion for /(x)=x. 


E. Integrals of Series Expansions 


If a function is represented by Eq. 8, its derivative with respect to x is not explicitly represented by the 
derivative of its component parts in the summation. The basis functions g n {x) have 0 slope on the interior 
of each segment and infinite slope at segment boundaries. For this reason, in dealing with any equation 
of mathematical physics involving derivatives of various orders, the derivative of highest order may be 
represented by Eq. 8 and subsequent lower-order derivatives of the function are obtained by integration. 
Thus, 


w~(x) = fx(x) = ^2 A n g n (x ) 


dx 


and 


pX 00 00 px 

f(x ) = f{x a ) +1^2 A ng n {s)ds = f(x a ) + ^ A n g n (s)ds 

J Xo, n=1 n= l -Ixo. 


An integrating matrix x(n,m) is defined such that 

pX 00 

/ g n {s)ds = ( x b - x a ) ^2 x{n, m)g m {x) 

J X n 1 


(18) 


(19) 


( 20 ) 


m = 1 


on the interval x a < x < Xb ■ The factor ( Xb — x a ) is introduced to keep x(n,m) dimensionless. The 
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coefficients x( n i m ) are given by 


X(n,m) 


(xb - x a ) 


- 3 / 2 E' 


-i) 


(i+l) 


g n {s)di 


dx 


(21) 


where the Xi are defined in Eq. 11. The first 16 rows and columns of the integrating matrix are presented 
in Eq. 22. It is seen to be sparse; the occurrence of non-zero elements is easily defined so that integrating 
transforms need not be especially costly, even for large n. 


x( n i m) = 


1 
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0 


-1 

32 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


( 22 ) 


It is now possible to represent f(x) as 

oo oo oo 

f(x) = f(x a ) + (x b - x a ) EE A n x{n,m)g m (x ) = E (*r) 

m — 1 n—1 m— 1 

where 

A rf i — {x b Xa)x (m, 7l)A n T <5l m{Xb X a ) ^ f(Xa) 


(23) 

(24) 


F. Derivatives of Series Expansions 

A differentiation operation — computing the series representation of ^ from the series representation for 
/( x) — may now be derived from Eq. 24. Differentiation of a function represented by a series truncated to 
2 P terms is implemented as follows. Rewrite Eq. 24 using the transpose of the integrating matrix to obtain 

Am — ( X b Xa)x (m, Tl)A n T 5i m (Xb Xa) ^ f(x a ) (25) 

where repeated index n indicates summation. Recall that A m are the coefficients of the series expansion for 
f(x) in Eq. 23 and A n are the coefficients for the series expansion of f x (x) in Eq. 18. If one truncates the 
series representations at n = 2 P in Eq. 25, then the inverse of the transpose of the integrating matrix may be 
calculated. The truncated series coefficients for f x (x) are now derived from the truncated series coefficients 

of /Or). 


A n = ( x b -x a ) 1 [x T (m,n)] 1 A m - Sim{x b - x a ) 1/2 f(x a ) (26) 

A n = (x b - £ a ) -1 £> p (n,m) A m - S lm (x b - x a ) 1/2 f(x a ) (27) 

where summation over repeated index m is implied and D p = [x T (m, n)] _1 with 1 < m < 2 P and 1 < n < 2 P . 
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Note that the elements of the integrating matrix \ are not a function of the truncation level of the series 
but the elements of the differentiating matrix T> p are a function of the truncation level n = 2 P . For example, 


V 2 


0 0 0 -8 
0 0-8 0 
0 8 0 -16 
8 0 16 32 


V 3 


0 0 
0 0 
0 0 
0 0 
0 0 
0 0 
0 16 
16 0 
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0 16 0 

16 0 0 
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0 -32 0 

32 0 -64 

0 64 128 


(28) 
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0 0 0 0 

0 0 0 32 

0 0 32 0 

0 32 0 0 

32 0 0 0 
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0 0 
0 0 
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-32 0 

0 0 
0 0 
0 0 
0 0 
0 0 
0 0 
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0 0 
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128 0 


0 -32 
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0 0 

0 0 

0 0 

0 0 

0 0 
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-64 0 

0 0 

0 0 

0 -128 

-128 0 

0 -256 

256 512 


(29) 


The sum of elements in any column is 2 P+1 , the sum of elements in all but the last row is — 2 P+1 , and 
the sum of elements in the last row is 2 2p+2 — 2 P+1 . The only non-zero element of the first column of the 
differentiating matrix of any order p is in the last row. 



X 


(a) C = 0. 

Figure 12. Series representation of cos(aj) computed 
the interval -3 < x < 3 for p = 5 . 



the derivative of the series representation of sin(aj) on 


For smooth functions it is observed that the presence of f(x a ) serves to drive the value of A^v towards 
zero in Eq. 27. As an example of this effect, consider again the sample problem of computing the series 
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representation of cos(x) as the derivative of the series representation of sin(a;) using the differentiating matrix. 
However, in this case change the interval to —3 < x < 3. If the contribution of /( x a ) is replaced by an 
arbitrary constant C and C is set to zero, then the differentiated series representation yields the comparison 
shown in Fig. 12(a) with L\ = 1.5097 and A 2P = —3.698. If one chooses C such that A 2 p = 0, then the 
comparison in Fig. 12(b) is obtained with L\ = 0.00091067. If C = f(x a ) = sin(— 3) then A 2P = —0.0108 
and L\ = 0.004425. It should be noted that if the definition of L\ in Eq. 13 used the difference of the series 
representation relative to the function rather than the absolute value of the difference then the representations 
in Figs. 12(a and b) have identical norms equal to 0.000069322. 
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(c) p = 8, full approximation 


(d) p = 8, limited approximation 


Figure 13. Series representation of Eq. 30 and its derivative using Eq. 27. 


If the function f(x) has a discontinuity or it is inadequately resolved by its series representation A n g n (x), 
then the representation of its derivative f x (x ) by the series A n g n (x) through Eq. 27 is observed to oscillate 
around the exact solution. This behavior is illustrated in the left side of Fig. 13 where the function has 
discontinuities at x\ and x 2 . 
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In Fig. 13, the original function is shown in red, its series representation is shown as a thick black line, and 
the derivative of the series representation from Eq. 27 is shown by the thin black line. The derivative of the 
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original function equals zero everywhere except at the points of discontinuity where 


lim 

e-J-0 


/ Qi + e) - f(x 1 
2e 


e) 


— > 00 


and 


lim 

e-s-0 


f(x 2 + e) - f(x 2 
2 e 


0 


— >■ —00 


(31) 


In all cases, the series representation of the original function is good, consistent with the earlier example 
of the Heaviside function showing no overshoots or Gibb’s phenomena. However, the series representation 
of the derivative of the function exhibits oscillations of constant magnitude over each piecewise constant 
segment of the original function that extend all the way to the left and right boundaries. The mean value of 
the oscillations equals zero over these piecewise constant segments. A representation that averaged this result 
over a segment length equal to 2dx p provides an excellent result. This averaging is easily accommodated by 
truncating the highest frequency components of the series from 2 P_1 < n < 2 P . The right side of Fig. 13 
shows this limited series representation of the derivative. The result is implemented by replacing T> p in Eq. 
27 with T>f im where 


D Urn( n ’ m ) 



for 

n < 2 P ” 1 

for 

n > 2 p- 1 


(32) 


The "spike" over the discontinuity has width equal to dx p for the original formulation and is twice as wide 
in the limited formulation. As p increases the width of the spike decreases and its height grows. Limiting 
coarsens the resolving ability of the original series by a factor of 2. The limiting is applied uniformly over 
the domain and will affect both smooth (well resolved) and discontinuous (poorly resolved) functions. 


III. Truncation and Prolongation 

Truncation (filtering) refers to the reduction or elimination of high-frequency (small-wavelength) struc- 
tures due to noise or incipient instability from a function represented by a Walsh series. Prolongation 
(sharpening) refers to the insertion of high-frequency (small-wavelength) structure into a Walsh series to 
smoothly replace eliminated noise or to serve as the first step of a wave-component-sequencing algorithm 
(equivalent to a grid-sequencing algorithm in a finite-difference context) . Truncation is effective in reducing 
Gibbs phenomena in discontinuous solutions and in suppressing instabilities in high-order formulations of 
derivatives. 

Consider the baseline distribution qo(t n ) = cos( 27 t t n ) with 1 < n < N and t n = (n — 0.5) /N where 
N = 2 Pt . Use a random number generator to add noise to the distribution such that qf{t n ) = qo(tn) + 
0.1(rn(f ra ) — 0.5) where m{t n ) is a random number in the interval 0 < m(t n ) < 1. Let 


2 p T 

Qf(t) = ^2 Af n9n{t) (33) 

n—1 

Figure 14(a) shows the noisy, unfiltered distribution in green plotted over the smooth baseline distribution 
in black. It also shows the time derivative of the noisy, unfiltered solution in blue plotted over the time 
derivative of the baseline solution in red. A single level of filtering ( jif = 1 ) is implemented by setting 
the wave components in the highest group number (recall Eq. 7) to zero ( p = 10 in this example). This 
action effectively averages adjacent values of qf{t 2 i-i) and qf(t 2 i-\) across these time intervals defined by 
1 < i <' 2 10_1 , but with fewer operations. Two levels of filtering (n/ = 2) are implemented by setting the 
wave components in the two highest group numbers to zero {p = 10 and p = 9 in this example). In like 
manner, this action effectively averages four adjacent values of < 7 /(t 4 i_i), <Z / (^ 4 i— 2 ) ■. qf(ta- 3 ). and qf(t 4 ^) 
across these time intervals defined by 1 < i < 2 10 ^ 2 . Each successive level of filtering (nj) using function 
trun below effectively averages another factor of 2 intervals to form a larger, filtering interval. 

2 (p T -r lf ') 

trun (q f (t),n f )= ^ A f n g n (t) (34) 

n = 1 

The baseline average value of qo{t) is accurately approximated at the center of the filtering interval if 
one assumes ( 1 ) the variation of the baseline solution is nearly linear across the filtering interval; and ( 2 ) 
the filtering interval is sufficiently long that noise fluctuations average to zero. In Fig. 14(b) a single level 
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of filtering creates a filtering interval of length 2~ < 4° _1 ) that shows marginal reduction of fluctuating noise. 
In Fig. 14(c), five levels of filtering creates a filtering interval of length 2 - ( 10-5 ) that averages noise over 2 5 
adjacent segments. Note that trun(<7/(i„), 5) (in green) shows a 32 stair-step pattern that is flat over each 
set of adjacent 2 5 intervals. The baseline function qo(t n ) (in black) is nearly linear over each of these steps. 

The five-level filtering in Fig. 14(c) has discarded 2 10 — 2 10 ^ 5 wave components (all but 32 out of the 
original 1024). The time derivative of this piecewise constant distribution from Eq. 27 is zero as shown by 
the blue line in the figure. (Note that Eq. 27 is not equivalent to a conventional Taylor series representation 
of a derivative in which a divided difference across the stair riser will always yield a large, non-zero result.) 
One would like to recover a smooth variation across the original 1024 segments without having to resort to 
generating a curve fit through the remaining 32 segments. 




(a) qf , unfiltered 


(b) qf , filter level 1, with prolongation 




(c) qf , filter level 5, without prolongation 


(d) qf , filter level 5, with prolongation 


Figure 14. The use of filtering and prolongation to reduce random noise. 


The fractal nature of the Walsh basis functions 1 introduces repeated patterns in the wave component 
magnitudes that can be exploited for the purpose of prolongation — replacing the truncated, high-wave- 
number components of the noisy solution with components that yield a smooth solution over the filtering 
interval. The repeated patterns are suggested in Figs. 7(b), 8(b), 9(b), and 10(b). Knowing to expect 
patterns makes it easier to find them. To this end, plot the wave components using a mapping that retains 
group number identity as follows: 
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to = 2 P + 1 — n for 2 P 1 < n < 2 P (35) 

Equation 35 plots wave components for each group in reverse sequency order starting from to = 1. It 
enables the superposition of wave components to better reveal underlying patterns as demonstrated in Figs. 
15(a-c). (Information regarding the sign of a wave component is lost here in order to focus on component 
magnitude using logarithms.) The base 2 logarithm of the Walsh series components of three simple functions 
on the interval 0 < t < 1 are plotted as a function of group number using Eq. 35. Elements of each group 
are plotted with a different symbol and color. These plots are referred to as fractal fingerprints (FFP) of a 
function over an interval. A one-unit change in these figures between component values in adjacent groups 
indicates a factor of 2 difference in magnitude. If a function is singularity-free over an interval, then the 
magnitude of adjacent groups differ by approximately a factor of 2. This trend is highlighted in Figs. 15(d f) 
in which the FFP is plotted after dividing the wave components in group p by a factor 2 P . 
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Figure 15. The fractal fingerprint (FFP) of Walsh series wave components. 


Examples of the FFP for functions that include singularities in the interval 0 < t < 1 are presented in 
Fig. 16. Note that successive groups still show repeating patterns; however, the magnitude of successive 
groups cannot be described by a single factor of 0.5. If the singularity is moved outside of the domain (for 
example, q(t ) = 1 /(t + 0.1)), the FFP shows good recovery of the 0.5 scale factor between groups. 

The repeated patterns in the FFPs are used as a first approximation to increase the resolution of a Walsh 
series representation of a function to the next highest group level. It has already been noted that the left 
half of the highest group number is already well-approximated by simply setting it equal to half the value 
of the preceding group. It will next be observed that polynomials of degree p exhibit a repeating pattern 
in their FFP in which the right half of the highest group number can be predicted from the values in the 
left half. (These repeated patterns for polynomials are consistently observed but not yet proved.) These 
patterns are exploited to fully construct a good approximation to wave components in group p completely 
from wave components in group p — 1 — a bootstrap approach to prolongation for the purpose of wave 
component sequencing or smoothing the structure that had been truncated in a filtering operation. 
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Figure 16. The fractal fingerprint (FFP) of Walsh series wave components for functions with singularities in 
the domain. 


The algorithm to define elements of group p from a Walsh series that terminates at group p — 1 is 
illustrated in Fig. 17. In Fig. 17(a), the scaled FFP for q(t) = t 6 is used to discern repeated patterns. First, 
note that the scaled elements of group 6 (black squares) in the left half domain match the scaled elements 
of group 5 (red delta). Next, note that the distribution in the far right quarter of group 6 (black squares, 
25 < to < 32) matches the distribution in the second quarter of group 6 (9 < to < 16) but is displaced 
vertically by one unit (a factor of 2 smaller). Next, observe that the distribution in the first half of this 
latest group (25 < to < 28) is repeated at the front eighth of the second half (17 < to < 20) but is again 
displaced vertically by one unit (a factor of 2 smaller). This fractal algorithm of transferring half of the 
last defined pattern to the opposite end of the remaining undefined domain and reducing the magnitude by 
factor 2 continues until the entire group is defined. 

A schematic of this prolongation algorithm is illustrated in Fig. 17(b). The logical flow goes from the top 
to the bottom of the figure. Red blocks represent existing patterns. Blue blocks represent target locations 
where the pattern in red is transferred with magnitude reduced by a factor f, in step i of the algorithm. 
Gray blocks represent locations in the new group that have been defined previously. White areas have yet 
to be filled in. Each step takes half of the latest defined subdomain and transfers it to the opposite end of 
the undefined domain, reducing magnitude by a factor /,; until the entire group is defined. The reduction 
factor fi is currently fixed at 0.5 based on the polynomial template function. Other values may be derived 
for other template functions, but such an option is beyond the scope of this paper. The current algorithm 
provides good approximations to the prolongation of singularity-free functions across a domain; errors tend 
to build on the last reconstruction steps, but these coefficients tend to be orders of magnitude smaller than 
other coefficients in the Walsh series. 

Figure 17(c) confirms that the reconstructed group 6 coefficients (red circles) obtained by prolongation 
from the group 5 coefficients match the exact group 6 coefficients (black squares) obtained from a Fast Walsh 


17 of 44 


American Institute of Aeronautics and Astronautics 


Transform of q(t) = t 6 . The sequential application of this prolongation algorithm from group 5 to group 10 
was used to recover a smooth baseline solution in Fig. 14(d) (green curve) from the truncated solution in 
Fig. 14(c). The time derivative of this reconstructed solution (blue curve in Fig. 14(d)) shows some small 
discontinuities associated with approximate nature of the bootstrap algorithm for general functions. 




(a) Scaled FFP for q(t) = x 6 including (b) Algorithm schematic to bootstrap (c) Comparison of reconstructed elements 
groups 0 < p < 6. The elements in group group p using only elements of group to original, untruncated elements of group 
6, black squares, are truncated from this p — 1. 6. 

series. 


Figure 17. Overview of prolongation algorithm applied to Walsh series representation of q(t) = x 6 starting 
with 6 groups. 


IV. Application Examples 

Three test cases provide examples on the use of the Walsh functions to solve partial differential equations. 
All of the sample problems are time-dependent. Two are non-linear. One involves a system of equations. 
All solutions utilize the Fast Walsh Transform (FWT) documented in Appendix A to produce the discrete 
solution values shown in the figures from the Walsh function wave components. 

The formulations in Section A are completely executed in Walsh wave component space. Boundary 
conditions are implicitly included in the formulation of derivatives using Eq. 27. The advantage of this 
approach is a robust, quadratic convergent algorithm. The Jacobians of every wave component with respect 
to changes in every other wave component and with respect to boundary conditions are automatically 
implemented using operator overloading and the chain rule. The disadvantage of this approach is that 
Jacobian matrices are not sparse and require large amounts of memory. Furthermore, two-point boundary 
value problems require the introduction of a second-derivative dissipation in order to accommodate two 
boundary conditions through Eq. 27. 

The formulations in Section C replace derivative formulations from Eq. 27 with traditional, finite- 
difference formulations derived from Taylor series. Finite-difference operations are enabled by a Walsh 
function shift in Eq. 36. 


shift(.F, inc, F bc ) 


' FWT (F) f 
FWT(F bc ) -> f bc 

fi ^ fi+inc 

if inc = 1 then f bc — > fi 

if inc = — 1 then f bc — > fx i 

^FWT(f) ->• shift 


(36) 


This equation uses (1) an FWT from wave components F a to discrete components /); (2) a shift of the discrete 
indices by ±1; (3) insertion of the appropriate boundary condition at the first or last index, respectively; 
and (4) an FWT back to wave components F a from discrete components fi. This approach enables explicit 
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inclusion of boundary conditions and does not require storage of the full Jacobian. More relaxation steps 
are required per time step than with the implicit formulation but convergence is orders of magnitude faster 
than the implicit formulation for problems with more than 1024 degrees of freedom. 


A. Test Case Formulation: Implicit Method 

1. Advection 

The linear, first-order advection equation is 


du 

dt 


du 

dx 


= 0 


(37) 


where c = 1 is the wave speed and u(x, 0) = uq(x) is the initial condition on the domain 0 < x < 1. The 
initial profile moves to the right. A periodic boundary condition is applied so that as the profile exits the 
right boundary it reemerges from the left. The test is designed to provide an easily evaluated metric to 
measure how well the initial profile is preserved. The initial profile in the demonstration case (see Fig. 18) 
is a sawtooth defined by 



f 0 

for 

0 < x < \ 


u 0 {x) = i 


for 

i < x < ^ 

(38) 


{ 1 

for 

! <*< i 




Figure 18. The initial profile for the advection test case, u$(x). 


The basis function series representation of u(x , t) is 

N a N t 

u(x,t) = ^2'52u(i,j)g i (t)g j (x) (39) 

i — 1 *= 1 

The indices i and r indicate a basis function component in time t. The indices j and a indicate a basis 
function component in space x. The coefficient U(i,j ) should never be interpreted as the value of u{x,t) at a 
specific point in the space-time domain. Rather, it indicates the magnitude of a basis function component of 
the global solution for u(x, t ) across the entire space-time domain. The number of basis function components 
considered in the x domain is N a = 2 Pa . The number of basis function components considered in the t domain 
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is N t = 2 Pt . It is important to include every member of a family p a and p T in non-linear problems to capture 
every basis function component that results from a multiplication of two or more functions. 

The basis function series representation of ^g^ = u x (x,t) is 

N a 

u x (x,t) =Y2^Z u x (i,j)g i (t)g j (x) (40) 

3 = i i = i 

The coefficients U x (i,j) may be expressed as a function of the coefficients U(i,j) using the integrating matrix 
from Eq. 24. 

N a 

U(i,a) = C x ^ U x (i, j)Xp a (j, a) + S(a, l)U x p(i) (41) 

3=1 

The factor C x is the length of the x domain. It is convenient in the case of a two-point boundary value 
problem to replace the last term of Eq. 24 with a constant of integration that will be solved implicitly with 
the other expansion coefficients. The constant of integration U x p(i) is most generally a function of time with 
basis function component i. The function 5 is Kronecker delta. Einstein summation notation is adopted to 
derive U x (i,j) as an explicit function of U(i,j) and U x p(i). 

U x (i,j) = C- 1 [D Pa (j, a)U (i, a) - V Pa (j, 1 )U x0 {i)} (42) 

Here T> Pa (j,a) equals the inverse of xj ct (a,j). 

The basis function series representation of au< ^^ = ut{x,t) is 

N a Nt 

Ut{x,t) ='^2'^2u t (i,j)g i (t)g j (x) (43) 

3 = 1 »= i 

The coefficients Ut(i,j ) are expressed as a function of the coefficients U(i,j ) and constant of integration 
Uto(j ) using the same algorithm defined in Eqs. 41-42. Therefore, 

Ut(i,j) = A" 1 P P Ahr)U(r,j) - V PT (h l)Ao(j)] (44) 

where the factor C t is the length of the t domain (the time step). A simple, first-order, backward difference 
approximation may also be expressed for the case p T = 0 and i = 1. 

Ut(l,j) = A” 1 P(l,j) ~ U t0 (j)} (45) 

A residual equation representation of Eq. 37 for the Walsh function component magnitudes is 

= U t {i,j) +cU x (i,j) (46) 

The amplitude of wave components the dependent variables, are solved such that each wave compo- 

nent R(i,j) goes to zero. A FORTRAN module, WALSH_TOOLS, 2 simplifies this solution. 

Code for Eq. 37 using WALSH_TOOLS is written 

if(N_tau > l)then 

resw(l) = intt (qlw(m) , f a=qlOw(m) , diff = . true . ) & 

+ wave_speed*intx(qlw(m) , f a=qlaw(m) , dif f = . true . ) 

else 

resw(l) = (qlw(m) - qlOw(m))/dt & 

+ wave_speed*intx(qlw(m) , f a=qlaw(m) , dif f = .true . ) 

end if 

All variable names in this example ending in letter w are of type walsh and include wave components 
and their Jacobians with respect to the wave components of the dependent variable. The function call 
intx(qlw(m) ,fa=qlaw(m) ,diff=.true. ) computes U x (i,j) in Eq. 42. In like manner, the function call 
intt (qlw(m) , fa=qlOw(m) , dif f= .true .) computes in Eq. 44. For the case with p T = 0 (N_tau = 1) 

the expression (qlw(m) - qlOw(m) )/dt computes Ut.(i,j) in Eq. 45. The functions intt and intx implicitly 
include the differentiating matrices T> Pt (i, r) and T> Pa (j, a), respectively. Note that for this linear, first-order 
equation, only a single Walsh function dependent variable qlaw(m) is required to satisfy a boundary condition 
in space, and qlOw(m) is required to satisfy an initial condition in time. Also note that index m refers to a 
subdomain. The solutions that follow use only one domain. 
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2. Burgers Equation 

The non-linear, second-order Burgers equation is 


du 1 du 2 d 2 u 
dt + 2 l)x = 


(47) 


where the diffusivity v > 0. The initial condition is a linear function uq(x) = —x on the domain — 1 < x < 1 
with boundary conditions u a (t) = u(— l,t) = 1 and Ub(t) = u(l,t) = — 1. 

Given these initial and boundary conditions, the solution profile evolves according to the local value of 
u. At any given time, points with a positive value of u move to the right and points with a negative value 
of u move to the left with speed u. This motion is dissipated as a function of diffusivity v. Large values of 
v evolve profiles that are nearly linear between the boundary values. Small values of v evolve profiles with 
an abrupt, shock-like transition from 1 to -1 at x = 0. Examples will be provided in Sec. B2. 

If the diffusivity is specified 0, then an artificial diffusion term must be included to maintain a stable com- 
putation while accommodating boundary conditions from both the left and right. The artificial dissipation 
is defined 


= 7){ dx p) 2 


du 

dx 


(48) 


where dx p is the smallest segment size used in the Walsh series representation of the solution. The use of a 
second-order, artificial dissipation here and in the Riemann problem that follows presents a simple solution 
for communicating information from opposite boundaries. Such communication was not an issue in the 
linear advection test case in which all waves travel from left to right. It is thought that a characterstic-based 
formulation of this problem could offer better accuracy as v — > 0, and the transition from u = 1 to a = -1 
occurs over a length smaller than dx p . For now, only the artificial dissipation is used to address this issue. 

Code for Eq. 47 using WALSH_TOOLS 2 is written 


if(N_tau > l)then 

resw(l) = intt (qlw(m) , f a=qlOw(m) , diff = . true . ) k 

+ intx(0.5_dp*qlw(m)**2-tauxw,fa=qlbw(m) ,diff=.true. ) 

else 

resw(l) = (qlw(m) - qlOw(m))/dt k 

+ intx(0.5_dp*qlw(m)**2-tauxw,fa=qlbw(m) ,diff=.true. ) 

end if 


Representation of the time derivative is exactly the same as discussed above for the advection equation. The 
shear term, tauxw = vdu/dx, was computed just above this block as 

uxw = intx(qlw(m) , fa=qlaw(m) , diff = .true . ) 
if (nu==0 . _dp) then 

dx2 = 0 . 5_dp*dx_ef f **2 
tauxw = dx2*absw(uxw) *uxw 
else 

tauxw = nu*uxw 
end if 


Note that if v = 0, then the artificial dissipation is added as defined in Eq. 48. Also note that this second- 
order equation has two boundary conditions that must be engaged. The Walsh function dependent variables 
qlaw(m) used in the definition of uxw = du/dx and qlbw(m) used in the definition of 3(0. 5u 2 — udu/dx)/dx, 
provide additional degrees of freedom required to satisfy the boundary conditions. The Walsh function 
dependent variable qlOw(m) used in the definition of du/dt provides the additional degrees of freedom to 
satisfy the initial condition. 


3. Riemann Problem 

The Riemann problem in one dimension defines the compressible gas dynamic flow that follows the breaking 
of a virtual diaphragm separating two constant initial states. The solution typically involves the propagation 
of a shock, a contact discontinuity, and an expansion fan. A particular instance of the Riemann problem 
was defined by Sod' and is frequently used in the validation of computational fluid dynamics codes. 
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The one-dimensional, time-dependent, compressible gas dynamic equations are 


dp dpu 
dt dx 
dpu d{p + pu 2 ) 
dt dx 

dpE dpuH 
dt dx 

V 

e 

H 

The need for artificial dissipation was discussed in t 


0 

(49) 

0 

(50) 

0 

(51) 

(7 - 1 )pe 

(52) 

E -i 

(53) 

E + P - 

(54) 


P 


previous section and is added in Eqs. 55 - 57. 


dpu 
~dt + 
dpE 

~df 


dp 

dt 

d_ 

dx 


d_ 

dx 


pu — Vi 


p + pu — 


d_ 

dx 


puH — v 3 


dp 

dx 

dpu 

dx 

dpE 

dx 


0 

0 

0 


The artificial diffusion coefficients defined in Eqs. 58 - 60 are of order ( dx p ) 2 . 




^2 


^3 


(< dx p ) 2 
( dx p ) 2 
{dx p ) 2 


d -l 

dx 

dpu 

dx 

dpE 

dx 


i d/Xp j ” 


( dx p ) 2 
- (dx p ) 2 


(55) 

(56) 

(57) 

(58) 

(59) 

(60) 


The equations are solved on the domain — 1 < x < 1 with 7 = 1.4. Constant initial conditions on the left 
x < 0) and right (x > 0) for the Sod test case' are given by 


Pl — 1 
Pl = 1 

u L = 0 


Boundary conditions are given by 


dp 

dx 


(-M) 


= 0 


pu(— 1, t) = 0 


dpE 

dx 


(-M) 


= 0 


PR = 0.1 
PR = 0.125 
Ur = 0 


dp 


(M) = 0 


dx' 
pu( 1, t) = 0 
dpE 


dx 


(M)= 0 


Code for Eqs. 55 - 57 using WALSH_TOOLS 2 is written 

if(N_tau > l)then 

resw(l) = intt (qlw(m) , f a=qlOw(m) , diff = . true . ) & 

+ intx(q2w(m)-qlxw,fa=qlaw(m) ,diff=.true. ) 
resw(2) = intt (q2w(m) , fa=q20w(m) , diff =. true . ) & 

+ intx(pw + rhou2w - q2xw,fa=q2aw(m) ,diff=.true. ) 
resw(3) = intt (q3w(m) , fa=q30w(m) , diff =. true . ) & 
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else 

resw(l) 

resw(2) 

resw(3) 

end if 


+ intx(q2w(m) *htw - q3xw, f a=q3aw(m) , dif f = . true . ) 

= (qlw(m) - qlOw(m))/dt 

+ intx(q2w(m) -qlxw, f a=qlaw(m) , dif f = . true . ) 

= (q2w(m) - q20w(m))/dt 

+ intx(pw + rhou2w - q2xw , fa=q2aw(m) , dif f= . true . ) 
= (q3w(m) - q30w(m))/dt 

+ intx(q2w(m) *htw - q3xw, fa=q3aw(m) , dif f= .true . ) 


& 

& 

& 


The main difference between this code sample and previous ones in this section is that there are three 
dependent variables and three corresponding residual equations to be solved. The reciprocal of density is 
required to compute velocity and pressure from the conserved variables. The reciprocal evaluation uses the 
Fast Walsh Reciprocal documented in Sec. VIII (Appendix B). 


B. Test Case Solutions with Implicit Method 

1. Advection 

The Walsh function solution of the advected profile after one cycle is compared to the exact solution in Fig. 
19(a). If a Courant number is defined as the ratio of the distance traveled by a wave in time step dt divided 
by the smallest segment size dx p , then the Courant number in this case is cdt /dx p = 2.55. The error norm 
after one cycle equals 1.66 10” 3 . Some oscillation is evident at the front foot of the profile. If the case is 
rerun through Eq. 34 with n.f = 1 to smooth oscillations, the error norm after one cycle increases to 1.37 
1CT 2 , and the profile is shown in Fig. 19(b). Truncation of wave component contributions for n > 2 i> °~ 1 
eliminates the oscillations at the expense of smoothing out the discontinuities at the tip and feet of the 
profile. 

If temporal degrees of freedom are exchanged for spatial degrees of freedom by decreasing p T to 0 and 
increasing p a to 10, then the Courant number increases to 10.23, and the profile appears in Fig. 19(c). 

In stark contrast to the previous case, if spatial resolution is sacrificed for increased temporal revolution 
by setting p a = 6 and p T = 6 and increasing the time step dt by a factor of 100, then a complete cycle 
is computed in a single time step, as shown in Fig. 19(d) with zero error norm for the linear problem. 
Zero error means that the advected profile exactly crosses the segment midpoints. The Courant number in 
this case is 63, using the previous definition. A more representative definition of the Courant number is to 
consider the distance traveled by a wave in the smallest temporal segment size c{dt) p divided by the smallest 
segment size in space ( dx ) p , in which case the Courant number is equal to 1. This result exhibits a resonance 
in that the computed profile exactly repeats every time step, so that the llnorm = 0 in a single relaxation 
step because the initial condition at the beginning of the time step equals the converged condition at the 
end of the time step. 


2. Burgers Equation 

A steady solution is obtained for t > 6.6 based on attaining an llnorm < 10 _1 ° in a single relaxation step 
following an advance in time. The steady solution at t = 10 is shown in Fig. 20(a). The segmented nature of 
the underlying Walsh function support is evident in the stair stepping appearance of the solution. The stair 
stepping is less evident in Fig. 20(b) in which p a = 8 and segment size is a factor of 4 smaller. The error 
norm for this problem is recorded in Table 2. With each increment in p a , the number of segments is doubled 
and the error norm decreases by a factor of approximately 4, indicating second-order accuracy relative to 
the smallest segment size. 

As the shock thickness decreases below the width of (dx) p with decreasing is, oscillations may be encoun- 
tered in the solution. These oscillations are eliminated by truncating the highest family of Walsh function 
contributions to the solution at the conclusion of a time step. An example is presented in Fig. 21, where 
the solution without (a) and with (b) truncation is shown when is = 0.001 and p a = 8. The error norm 
without truncation is 7.27 10^ 3 and with truncation is 1.94 10 3 , indicating roughly a factor of 4 decrease 
in error. It is a subtle but important point to note that the highest order Walsh functions still contribute 
to the lower-order non-linear solution retained after truncation through the self-mapping property under 
multiplication. 
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(a) p a = 8, p T = 2, truncate = 0, dt = 0.01 



(c) p a = 10, p T = 0, truncate = 0, dt = 0.01 

Figure 19. Advection test problem showing profile 
function series solution. Red is the exact solution. 



(b) p a = 8, p T = 2, truncate = 1, dt = 0.01 



(d) Pa =6 , p T = 6, truncate = 0, dt = 1. 

one complete cycle at t = 1. Black is the Walsh 


Table 2. Error norms for Burgers equation with v = 0.1. 


Pa 

error norm 

error norm ratio 

3 

1.28 1CT 1 

- 

4 

1.33 10~ 2 

9.62 

5 

2.60 10~ 3 

5.12 

6 

5.89 10" 4 

4.41 

7 

1.40 10~ 4 

4.21 

8 

3.45 10~ 5 

4.06 

9 

8.92 10~ 6 

3.87 
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Figure 20. Effect of p a in case of v = 0.1, p T = 0, and p domain = 0 for Burgers equation. Black is the Walsh 
function series solution. Red is the exact solution. 


0.5 


3 0 - 


- 0.5 - 


Exact 

Walsh Function 


- 1 - 

-1 - 0.5 


0.5 


3 0 



- 0.5 - 


Exact 

Walsh Function 


0.5 


(a) truncate = 0 


(b) truncate = 1 


Figure 21. Effect of truncate in case of u = 0.001, p a = 8, p T = 0, and Pdomain = 0 for Burgers equation. Black 
is the Walsh function series solution. Red is the exact solution. 


3. Riemann Problem 

The simulation is run with a time step dt = 0.002 with eight sub-domains. The Walsh series representation 
of each dependent variable includes 128 terms (p Q = 7) in each sub-domain. Each sub-domain overlaps its 
neighbor by dx p (Eq. 2). The solution for density at t = 0.42 is presented in Fig. 22. The expansion moving 
to the left is evident for — .5 < x < 0. The contact discontinuity at x ~ 0.35 and the shock at x ss 0.7 move 
to the right. The constant states between the shock and the contact discontinuity, and between the contact 
discontinuity and the head of the expansion, are in excellent agreement with the exact solution. The tail of 
the expansion at x ~ —0.5 is slightly rounded. The contact discontinuity appears more dissipated than the 
shock. The largest differences are thought to be associated with the direction of characteristics approaching 
the discontinuities. 
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Figure 22. Density profile at t = 0.42 for Riemann problem. Black is the Walsh function series solution. Red 
is the exact solution. 


C. Test Case Formulation: Explicit, Upwind Method 

The essential difference between the explicit formulations in this section and the implicit formulation defined 
previously in Sec. A is use of the shift function (Eq. 36) to enable explicit treatment of boundary conditions 
in the formulation of spatial and temporal derivatives - replacing Eq. 27. Additional advantages include 
the ability to introduce accurate formulations higher than second-order and eliminating the requirement for 
explicit artificial dissipation for inviscid problems. In essence, numerical tests indicate that the numerical 
dissipation implicit in the upwind formulation is sufficient to provide good shock capturing. 

A generic formulation of the test problems from the previous section follows. Only one space dimension 
is considered, but extension to multiple space dimensions follows identical patterns. 



dq df 

dt dx 

= 0 

(61) 


dq .dq 

dt dx 

= 0 

(62) 


d q , RAR i5q 

Sf +RAR ^ 

= 0 

(63) 

| + iR(A + |A|)R->g 

+ iR(A-|A|)R->g 

= 0 

(64) 


dq + dq 5q 

Ot dx + dx 

= 0 

(65) 


dq Sf + di~ 

dt dx dx 

= 0 

(66) 
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Here A is the Jacobian of f with respect to q, R is a right eigenvector matrix of A and A is the diagonal 
eigenvalue matrix of A. Note that formulations based on Eq. 65 showed poor resolution of jump conditions 
across discontinuities in approximations higher than lst-order whereas Eq. 66 was required for conservation. 
The differentials ^ across the interface between segments i and i + 1 are computed 

^1/2 = ±A± (qi=Fi/2) (q» - q*=Fi) ( 67 ) 

Sf w,-i /2 = A+ (qt« - 1 / 2 ) (q w - shift (q w , l,q w ,6 c ,o)) (68) 

5f“ 1/2 = - A_ (q w, 1 / 2 ) (q^ - shift(q w ,-l,q Wi6CiJV +i)) (69) 

The differential expressions in Eqs. 68 and 69 are written with subscript w to emphasize that the underlying 
Walsh series representation is equivalent to the discrete formulation in Eq. 67. The interface average, q, is 
defined such that Roe’s Property U 8 is satisfied, in particular, 


A (qu+i/ 2 ) (q*+i - q») = (f»+i - f i) 


(70) 


The derivative will be approximated by a backward difference formula and will be approximated 
by a forward difference formula to inject an accurate zone of dependence and to provide strong diagonal 
dominance for the explicit, inviscid solver as follows: 


<9f±\ 

Ih ) 


i 


/df± 

\ dx 


l 




|'(3f±-4f± 1 +f± 2 )/(2d I ) 

< (llff - 18^! + 9f± 2 - 2f^g )/(6dx) 

{ (25f± - 48f± x + 36f± 2 - 16f± 3 + 3f± 4 )/(12 dx) 

[ (3^i/2 -5f± 3/2 )/(2dx) 

| ( 11<5f ^l/2 - 7(5f ^3/2 + 2 ^ 5/2 )/(6dx) 
i(25<5f±i/ 2 - 23<5f± 3/2 + 135f* 5/2 - 3fl* 7/2 )/(12da:) 

((6f± 1/2 -6f± 3 / 2 )/(2d X ) 

| ( 5<5f ^l/2 - 75f ^3/2 + 2 ^ 5/2 )/{6dx) 

i(13^ Tl/ 2 - 2351*3/2 + 13^5/2 - 35f^ 7/2 )/(12 dx) 


+ 0(dx 2 ) 

+ 0(dx 3 ) (71) 
+ 0(dx 4 ) 

+ 0(dx 2 ) 

+ 0(dx 3 ) (72) 
+ 0(dx 4 ) 

+ 0(dx 2 ) 

+ 0{dx 3 ) (73) 
+ 0(dx 4 ) 


[(< T l/2-< T 3/ 2 )/( 2 ^) +0(dx 2 ) 

(5< t i/2 - 7 < 2 F 3/2 + 2 < t5/2 )/(6^) + 0(dx 3 ) (74) 

l(13<5f^ Tl/2 -23(5f± T 3 /2 + 13Jf± T5/2 -3Jf^ 7/2 )/(12dx) + 0{dx 4 ) 


where dx is the smallest segment size in the Walsh series representation defined by Eq. 2. The shift function 
(Eq. 36) is used to compute all elements of Eq. 74 and explicitly introduces boundary conditions from the 
left for Jf + and from the right for 8i~ . Again, the subscript w emphasizes that these terms are represented 
by the Walsh series that span the entire domain (all i). The terms to the right of the brace in Eq. 74 are 
treated separately as second-, third-, and fourth-order corrections (dfj^, o = 2, 3, 4, respectively) to the 
baseline, first-order derivative formulation. The magnitude of these corrections tend to be small except near 
discontinuities in W where they exhibit a relatively large magnitude oscillation. This behavior is clipped 
locally in discrete, physical space through clip(Jf 1 ± 0 ). 


dip(5f±) = < max[0, 1 - ( 1 
[ max[0, 1 — (— 1 



for 

frefYWto 

for 

5f to/freffWt 

for 


|5ftl< f ref 

6t to > fref ( 75 ) 

< -fref 


where f re f = f %, / XAI /i , o I / A * has been used herein and /° e j is a user defined clipping factor relative to 
the average value over the domain. 

Finally, the third- and fourth- order corrections to spatial derivatives are filtered to suppress instabilities 
using Eq. 34. The temporal derivative dq/dt is computed directly from Eq. 74 by substituting <5q„-i/2 for 
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<5f+ 1/2- It is not required to process the high-order corrections for temporal derivatives with the clipping 
function clip or the filtering function trun. 

The Walsh series representation of q w is obtained by driving the L\ norm of the residual of Eq. 66 to 
less than 10~ 10 over sub-iteration m for each time step as follows: 


= df^_ 94 

dt dx dx 


(76) 


m+l _ m 

1 w 1 w 


dr w 

dq w 


(q 


m+l 

w 


q») 


(77) 


(78) 

where the Jacobian of the Walsh series representation of the residual with respect to the Walsh series 
representation of the dependent variable is defined by 


m+l = m _ 

T.19 ^±W 


dr u 

9q u 


~9r w ] 9 / 9q,„ \ 9 ( 9f ,+ \ 9 / df w \ 

.AhJ dq w V dt ) 8q w { dx ) dq w { dx J 

and each component of the Jacobian is approximated by 


A / 

{ dq w 

dq,,, ' 

v dt 

A , 

AC 

5q w 

l dx 

A | 

( dC 

dq w 

^ dx 


c r 


dt 



d 

dq w 


- C “ A + 

d 

dq w 

(■ K,+l/J dx ) f 

~ Ca a 

dx ”’>+ 1 /2 


(79) 


(80) 

(81) 

(82) 


The relaxation coefficients c T and c a are usually set equal to the leading coefficient of <5++ / 2 /(dx) in Eq. 72. 
The Jacobian in Eq. 79 is a (3x3) for the one-dimensional Riemann problem where each element includes 2 P 
components of the Walsh series defining the element. Its inverse is calculated using Cramer’s rule for this 
small system. In contrast, the linear system in the implicit formulation has size (3 (2 P ) x 3 (2 P )), neglecting 
implicit treatment of boundary conditions. 


D. Test Case Solutions with Explicit Method 

1. Advection 

Solutions for the linear, first-order advection equation (Eq. 37) are repeated here with the explicit method. 
The wave speed c is a constant equal to 1; consequently A + = 1 and A~ = 0. In addition to the sawtooth 
initial condition (Eq. 38) a singularity-free Gaussian initial condition (Eq. 83) is tested to confirm the order 
of accuracy. 

uq(x) = exp(— 10(x — 0.5) 2 ) (83) 

Figure 23 summarizes results for the advection of a Gaussian profile. All cases examine the distortion of 
the initial profile after 20 cycles (t = 20) through the domain. A periodic boundary is specified so that the 
profile exiting the right boundary reemerges through the left boundary. None of the simulations required 
any clipping operations (/° e j = 10000) because the profile is continuous. The ratio of time step to minimum 
segment size is fixed at 0.9216. The domain is only one segment in time (p T = 0) in these tests. 

A progression of solutions with increasing resolution in space and time is shown in Figs. 23(a-d) for the 
third-order formulation with a single filtering step. An oscillation at the tail of the profile grows large after 
20 cycles for the largest segment size (p = 6). A single doubling of the number of segments (and halving 
the time step) still shows a slight oscillation at the tail of the profile after 20 cycles for p = 7 in Fig. 23(b). 
The oscillation is eliminated in subsequent refinements in Figs. 23(c-d). The second-order formulation with 
p = 9 did not require any filtering steps in Fig. 23(e); however, preservation of the initial profile is not 
quite as good as the corresponding third-order formulation. Solution accuracy is quantified in Fig. 23(f). 
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Third-order accuracy is confirmed with a slope 3 reduction in the log of the error norm versus the log of 
the segment size. Second-order accuracy is confirmed with a slope 2 reduction in the log of the error norm 
versus the log of the segment size. Fourth-order simulations at these conditions required a minimum of two 
filtering steps to suppress instabilities at the tail of the profile over large times. That filtering precluded a 
slope 4 reduction in error norm. In fact, the third-order result with a single filtering step was consistently 
more accurate than the fourth-order formulation with two or more filtering steps. 





(a) p = 6, third-order, rif = 1 


(b) p = 7, third-order, rif = 1 


(c) p = 8, third-order, rif = 1 



Figure 23. Advected Gaussian profile after 20 cycles through the domain. 


Advection of the sawtooth profile is reexamined in Fig. 24 using the explicit method. Second- and 
third-order formulations are presented in Figs. 24(a and b) with p = 9. Both cases use a CFL (dt / dx m i n ) = 
0.9216. The third-order formulation requires two filtering steps to suppress oscillations. The CFL is raised 
to 2.048 in Fig. 24(c) for the third-order formulation, but an extra filtering step was required to suppress 
oscillations. In Fig. 24(d), it is confirmed that clipping (/ 9 e j) is an ineffective substitute for filtering to 
suppress instabilities for this Co continuous profile. Though not presented here, additional tests of both 
second- and third-order formulations only produce slightly better than slope 1 (first-order) error reduction 
in this problem that is C\ discontinuous. 

2. Burgers Equation with Source Term 

The method of manufactured solutions 9 is used to modify Eq. 47 and verify order-of-accuracy for an unsteady, 
non-linear problem. A simple change to the initial condition enables an interesting unsteady problem with 
a moving shock with changing direction. Because the explicit method does not require a second derivative 
dissipation to close the problem (as with the implicit method), only the inviscid form (v = 0) in Eq. 47 is 
considered. 

The manufactured solution for this test is 

u(x, t) = 1 + sin(:r + 2 1) (84) 
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X 

(a) p = 9, second-order, rif = 0, f® e f = 10 4 , At = 0.0018 



X 


(b) p = 9, third-order, n f = 2, /° e/ = 10 4 , At = 0.0018 




(c) p = 9, third-order, rif = 3, /J? e ^ = 10 4 , At = 0.004 (d) p = 9, third-order, rif = 1, = 4, At = 0.0018 

Figure 24. Advected Delta profile after 20 cycles through the domain. 


Substitute Eq. 84 into Eq. 47 (with v = 0) to compute the requisite source term to obtain the manufactured 
solution. 


du 1 dn 2 , 

a +2^ = s(x ' 4) 

9(1 + sin(x + 2t) 1 d [1 + sin(:c + 2t)] 2 

dt 2 dx s x,t 

2 cos(a; + 21) + [1 + sin(x + 21)] cos(:r + 21) = s(x, 1) 


(85) 

( 86 ) 
(87) 


Two initial conditions are tested. In the first, u(x, 0) = 1 + sin(x) enables verification that the manufactured 
solution is recovered by the explicit formulation to the specified order of accuracy. In the second, u(x, 0) = —x 
produces a solution that approaches the manufactured solution through a shocked, intermediate state. 
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(a) L\ norm, second-order, exact initial (b) L ± norm, third-order, exact initial 
condition condition 


(c) p = 8, t = 2.00, u(x , 0) = —x 





(d) p = 10, t = 0.01, u(x , 0) = —x 


(e) p = 10, t = 0.50, u(x, 0) = — a; 


(f) p = 10, t = 0.75, u(x , 0) = —x 





(g) p = 10, t = 1.00, u(x, 0) = —x 


(h) p = 10, t = 1.50, u(x , 0) = —a; 


(i) p = 10, t = 2.00, u(a;,0) = —a: 


Figure 25. Burgers equation with source term from Method of Manufactured Solutions. 


In the explicit formulation 


UlAI, — 1/2 — 

I (q™ + shift (q,„, 1, q„,, bCi o)) 

(88) 

q^+1/2 = 

2 (q» + shift (q^,, -1, q w ,b c ,N+i)) 

(89) 

A ± = 

2 (q™,Tl/2 ± qu',Tl/2 ) 

(90) 


Verification of the convergence rate with the order of accuracy for the first initial condition yielding a 
continuous solution for all time is shown in Figs. 25 (a and b). The L\ norm at t = 2 is recorded with 
the second-order formulation showing a slope of 2 and the third-order formulation showing a slope of 4. 
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The third-order results include a single filtering sweep. Fourth-order formulations (not shown) required two 
filtering sweeps that prevented a slope 4 convergence rate. 

Some simulations with p a = 8 were executed with finer resolution in the time domain (p T = 1,2,3). The 
time step was increased by a factor of 2 Pt in order that the minimum time step used to compute a CFL was 
consistent with the other runs. Note that the error norms for all of the p a = 8 cases are consistent. 

Figures 25(c-i) present results for the second initial condition that yields a moving shock as the solution 
evolves in time. The red curve in these figures shows the manufactured solution and the black curve shows 
the computed solution at the same point in time. All of the computed solutions use the third-order accurate 
formulation with a single filtering step (n/ = 1) and a clipping parameter (/° e y = 4). Fig. 25(c) can be 
contrasted to Fig. 25(i) to see the sharper jump that is enabled by better spatial resolution. In all cases, 
there is usually only one segment and never more than two segments in the captured shock. The temporal 
evolution of the profile is evident in Figs. 25(d-i). The linear, initial profile is still evident at t = 0.01 in Fig. 
25(d). At t = 0.5 (Fig. 25(e)), the computed profile starts to collapse on to the reference profile at the left 
boundary. The profile is beginning to steepen at x = 0.6 because of left moving waves (u < 0) associated 
with the initial condition. At t = 0.75 (Fig. 25(f)), the discontinuity is being pushed to the right where the 
magnitude of u to the left of the jump exceeds the magnitude of u to its right. By t = 1.0 (Fig. 25(g)), the 
shock is just starting to move back to the left as the source term is pulling the magnitude of the left side of 
the jump condition to lower values. The shock continues to be pulled to the left in Figs. 25(h and i) as the 
reference solution approaches a minimum value before beginning to rise again. Though not shown here, the 
shock continues to the left until stopping at t = 3.08, x s=s —0.4, and then sweeping quickly back to the right. 
Other reversals of direction occur at t = 4.56 and t = 5.44 until the computed solution collapses completely 
on the manufactured solution for t > 6.2. It is observed that the clipping and filtering functions maintain a 
very sharp jump condition without oscillation or overshoots. 

3. Riemann Problem 

The explicit formulation of the Riemann problem follows for all quantities defined in Eqs. 66-69. All of the 
operations involve arguments defined by a Walsh series and use operator overloading in their implementation. 


[P,u,H] wSt+ 1 
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The Sod test case was discussed previously in Sec. B3, and its solution by the implicit formulation was 
presented in Fig. 22. Corresponding results with the explicit formulation are shown in Fig. 26. In this figure, 
the computed solution in black is compared to the exact solution in red. The simulations show results with 
9 < p < 12 and second- to fourth-order formulations of the derivatives with time step dt = 0.001. The L\ 
error norm for these cases is dominated by the solution at the moving shock and contact discontinuity and 
so error reduction as a function of grid refinement is first-order. In all cases, the shock (x = 0.74) is captured 
with fewer mesh points than the contact discontinuity ( x = 0.38). The tail of the expansion near x = 0 shows 
a sharp change in slope. The head of the expansion near x = —0.45 shows a rounded transition compared 
to the exact solution. It is difficult to discern any significant differences between second-, third-, and fourth- 
order accurate simulations where the only non-constant state is in the expansion region (—0.45 > x > 0). 
Increasing resolution has the most significant effect on solution accuracy as seen by comparing Fig. 26(f) to 
any of Figs. 26(a-e). Still, none of the simulations repeat the exceptionally sharp shock captures (only one 
interior segment) as achieved with Burgers equation in Figs. 25(f-i). Clipping the high-order correction is 
effective in eliminating overshoots in all cases. A single level of filtering the high-order corrections is required 
in the third- and fourth-order simulations to suppress instabilities. 





(a) p= 9, second-order, /(A = 0.5, nj = 0 (b) p = 9. third-order, /® e , = 0.5, nf = I (c) p = 9, fourth-order, //(^ = 0.5, n.j = I 



(d) p = 9, fourth-order, /° ,=0.5, np = 2 (e) p = 10, fourth-order, f® e j= 0.5, ny = 2 


(f) p=12, second-order, f® e f= 0.1, nf=0 


Figure 26. Riemann problem with Sod initial conditions at t = 0.4. 


The initial conditions across the diaphragm in Fig. 27 correspond to a pressure and density jump 
corresponding to a Mach number of 8.5. (Tests up to Mach 50 conditions have been generated with equivalent 
results.) These conditions produce a more challenging simulation in that a very narrow domain between the 
shock front (x « 0.3) and the contact discontinuity (x « 0.2) in Fig. 27(b) creates a slug of high density gas 
that is more than a factor of 2 denser than in the undisturbed state to its right or the post-expansion state 
to its left. The domain between the contact discontinuity and the tail of the expansion is slightly supersonic. 
The expansion from (—0.3 < x < 0) goes from subsonic to supersonic conditions. A small expansion shock 
is evident near x = 0 in Figs. 27(b-c). Numerical tests (not shown) indicate the jump is dissipated by 
setting | it ± c\ to be greater than c/4 in Eq. 64 - a simple implementation of the entropy fix. 10 Note that 
the simulation tracks the motion of all critical points (shock front, contact discontinuity, tail, and head of 
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the expansion) while running in explicit mode with a CFL = 2.13. 



(c) t = 1.0 (d) t = 1.5 

Figure 27. Riemann problem with p L = 1 , Pl = 6, pj j =0.01, and pp = 1 with At = 0.005 and A ma;c At/dx m i n = 
2.13. Numerical parameters are p a = 10, third-order, /® e y = 0.1, rif = 1. 


Table 3 compares the time to execute 10 time steps for the Riemann problem using the implicit and 
explicit formulations. Times are based on a compilation using a gfortran compiler with a -OO option on a 
2.93 GHz Intel Core 2 Duo processor. No attempt has been made to tune either the implicit or explicit path. 


Table 3. CPU time to advance 10 time steps. 

p a Implicit time Explicit time Iter /step, implicit Iter /step, explicit 


6 

12.933 

0.485 

2 

5.2 

7 

75.421 

0.733 

2.9 

5.2 

8 

377.712 

1.450 

3 

5.2 

9 

1691.462 

3.438 

3 

6.1 


The average number of relaxation steps per time step to drive the residual below 10 10 are presented in the 
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last two columns. The explicit method requires more relaxation steps per time step to achieve the target 
residual than the implicit method. However, the approximate Jacobian of the explicit method requires 
2 P fewer elements than the exact Jacobian of the implicit method thus enabling orders of magnitude of 
time-saving for p > 7. 


V. Summary 

Walsh functions form an orthonormal basis set consisting of square waves. The discontinuous nature 
of square waves make the system well suited for representing functions with discontinuities. A review of 
Walsh function fundamentals, starting with their derivation and including descriptions of algorithms to 
form their products, reciprocals, integrals, and derivatives is provided. Details of a Fast Walsh Transform 
in multi-dimensions is provided, which enables transforms from discrete variable space to wave space and 
back in order N log(iV) operations. A Fast Walsh Reciprocal is defined that enables the inversion of a Walsh 
Symmetric Matrix in order N log(-ZV) operations. This inversion is required in the evaluation of the reciprocal 
of a function represented by a Walsh series. 

Because Walsh functions can be derived using a fractal partitioning of a domain, 1 it is not surprising 
(perhaps even expected) that Walsh series representations of functions also exhibit recursive, self-similar 
patterns in their wave components. A remapping of Walsh function components is introduced here that 
makes it easier to identify these recursive patterns. The resulting figure is defined as a fractal fingerprint 
(FFP). The FFP of several functions that are singularity- free in a domain reveal a simple factor 2 scaling 
between components in adjacent groups. The FFP of functions that include singularities in the domain also 
show recursive patterns but the relation is more complex. These patterns are observed but not proved and 
more work is required to understand their potential utility in solving partial differential equations (PDEs). 

Three test problems focusing on unsteady simulations documented herein include linear advection, Burg- 
ers equation, and a Riemann problem. The first case provides examples of a profile with discontinuities in 
slope that are advected across the domain. The other two cases explore shock-capturing and the ability to 
accurately track a moving discontinuity. Two algorithms are developed to simulate these test problems. The 
first, a purely implicit algorithm, solves for the Walsh series components individually. While excellent con- 
vergence is achieved and the Jacobian of the solution is automatically generated through the use of derived 
types and operator overloading, the size of the Jacobian rapidly becomes intractable (on a single processor) 
for more than 2 13 degrees of freedom. Consequently, a second, explicit algorithm is derived in which the 
solution is focused on an update to the entire series representation of the function as a single entity - and 
not on individual wave components as if they were independent variables. Thus, a one-dimensional Riemann 
problem with 1024 degrees of freedom would require a linear solve of a (3072x3072) system where each 
element of the system is scalar in the implicit formulation. By contrast, the explicit formulation requires a 
linear solve of a (3x3) system where each element of the system includes 1024 Walsh series components that 
define the dependent variables over the entire domain. 

This paper presents the initial exploration of algorithms based on Walsh function series to solve un- 
steady, non-linear PDEs. The use of truncation to suppress instabilities and prolongation to serve as the 
wave-component counterpart of grid sequencing are demonstrated. Third-order solutions using the explicit 
formulation are verified using the method of manufactured solutions for a non-linear Burgers equation and 
comparing to the exact solution for advection of a Gaussian profile. Courant numbers (CFL) up to 2 have 
yielded accurate results for time dependent problems. Tests involving higher CFL numbers (greater than 5) 
are stable with filtering of high-order corrections but deform the solution profile, defeating the purpose of a 
time accurate simulation. A longer term goal to gather metrics comparing the evolving algorithms based on 
Walsh function series to more traditional approaches will follow from the work presented herein. 
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VII. Appendix A: Sequency-Ordered Fast Walsh Transform in 

Multi-dimensions 


Throughout this paper, Walsh functions are represented by g. Subscripts indicate a component of the 
Walsh series as a function of independent variables ( x , y , z, <). This representation is maintained to emphasize 
the origins of this algorithm as an orthonormal series representation of the solution. However, the mechanics 
of Fast Walsh Transforms (FWT) from discrete values of a function to its wave/sequency components and 
back makes explicit evaluation of g unnecessary provided that the number of components is a power of 
2 and the discrete functional values are uniformly distributed in computational space. The FWT enables 
transforms in O (N log(lV)) operations as opposed to 0(N 2 ) operations as would be required in a conventional 
evaluation of the series. 

Equation 102 presents the expansion of a function / of four independent variables (three space and one 
time) in terms of a Walsh function series . 1 In this equation, N a , Np, N 7 , and N T are the total number of 
Walsh functions used to represent the function in the x, y, z, and t directions, respectively. F a ^^ T are 
the wave number components of the expansion, where a, /?, 7 , and r are integers defining the number of 
segments in each Walsh function component direction. 

N t N 1 Np N a 

f(x,y,z,t) = EEEE -^ck,/3,7 ,r9a {x)gp {y)d^{z)g T (t) (102) 

t—1 7=1 /3—1 a—1 


Equation 103 defines the domain over which Eq. 102 is valid. This equation also sets up a uniform 
discretization of space based on the smallest segment size of the Walsh function in each respective direction. 
Here i, j, k, and n are integers defining the uniform discretization of space in the x, y, z, and t directions, 
respectively. 


Xo < x < XN a 
Vo<y< Vn p 
Zq<Z< z Ny 
to < t < t N T 


Ax = (x Na - x 0 )/N a 
Ay = ( y Nf> - yo)/Np 
A z = ( z Ni - z 0 )/N 7 
At = (: t Nr - to)/N r 


Xi = x o + iAx 
Vj = yo+jAy 
Zk — zq + kAz 
t n = to + nAt 


1 < i,a < N a 
1 <j,P<N 0 
1 < k, 7 < Noy 
1 < n, r < N t 


(103) 


Equation 104 defines the integral average of f(x, y, z, t) in the discretized element defined by Xi- \ < x < 
x i, Vj- i < y < y : j , Zk~i < z < Zk, and t n _i < t < t n . The representation of the discretized function fi.j t k,n 
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as a Walsh function series in Eq. 105 follows from Eq. 102. 


Ji,j,k,n 


1 


r z k rVj 


f(x,y,z,t) dx dy dz dt 


' Zk - 1 J Vj-i J Xi-i 


AxAyAzAt J tn l 

N t N 7 Np N a 

= EEEE ^a,/3,7 ,r<7a (,Xi— 1 / 2 ) yyiddj — 1/2)9^ (Zk—l/ 2 )yr (tn — 1 / 2 ) 

T— 1 7—1 /3=1 a=l 


(104) 

(105) 


It is convenient to think of terms like x^_ i / 2 as the point halfway between 27 _i and Xi . Because the 
Walsh function g a (x ) is piecewise constant over a segment, and all segment lengths are integer multiples (2 P ) 
of the smallest segment size, terms like Xi^i /2 may refer to any point between x,;_i and 27 . 

Equation 106 defines the Walsh function wave number components F a _p^_ T as a function of f{x,y,z,t). 
The definition of F a ,p,~t, T as a function of the discretized values fi,j,k,n in Eq. 107 follows from Eqs. 104 and 
106. 


r*JV T r z JV 7 I'V^I 3 r X N a 

Fa, 13 , 7 ,r = / / / / f(x, y, z, t)g a (x)gp(y)g 1 (z)g T (t)dxdydzdt (106) 

Jt 0 J z 0 Jyo J x 0 

N t N-, Np N a 

= EEEE fi,j,k,n9a(x i _ 1 /2)g0{yj-l/2)9 1 {z k _ 1 / 2 )gr{tn-l/2)^xAyAzAt (107) 

n— 1 k—1 j= 1 2—1 


If one has all of the elements F a j g i7jT - in Eq. 105, then any single element of fi,j,k,n can be computed in 
47V multiplications and N additions (assuming Walsh functions are pre-computed) where N = N a NpN 7 N T . 
Consequently, all N elements of fi,j,k,n can be computed in 0(N 2 ) operations. In like manner, if one has all 
N elements of fi,j,k,n in Eq. 107, then all N elements of F a ^ iJ:T can be computed in 0(N 2 ) operations. In 
the special case where each dimension is an integer power of 2 (N a = 2 Pa , Np = 2 Pp , N 1 = 2 P ~< , N T = 2 Pt ), 
a fast transform from all N of F a p ~ T to fij t k,n or back can be implemented in 0(pN ) operations where 
integer p = p a + pp + + p T - Closure under multiplication has this same requirement, 1 and so it adds no 

extra burden in the solution of differential equations. 

Note the similarities between Eq. 105 and its inverse transform in Eq. 107. If one executes a change in 
variable from (i,j, k , n) to (a, fi, 7 , r) in Eq. 107, then it has the same formulation except for the normalizing 
factor AxAyAzAt. and a transpose of subscripts in the Walsh function factors. Inspection of Fig. 7 of the 
preceding paper 1 on this topic shows that terms like g a (x.i_ i/ 2 ) are equal to their their transpose, 9 i{x a _ 1 / 2 )- 
Indeed, the matrix < 70 ( 27 - 1 / 2 ) is equal to its own inverse within a normalizing factor. 

Sequency-ordered (by increasing number of segments) Fast Walsh transforms are well documented for 
a function of a single variable . 11 ’ 12 A simple extension to multi-dimensions is defined here. A recursive 
algorithm 12 that does not require bit-reversal is modified here to accommodate multi-dimensional functions. 
This algorithm retains the advantage of producing its own inverse to within a normalizing factor. 

A schematic of the algorithm is presented in Fig. 28. The fast transform implemented here works 
sequentially through the x, y, z, and t dimensions. Scaled discrete functional values of fij,k,n (or of Walsh 
function component factors F a: p t 7jT ) are stored in a one-dimensional array f(t) where 


l 


i + N a (j — 1 + Np(k — 1 + N 7 (n — 1))) for transform from fi,j,k,n to F a> p tljT 
a + N a ((3 — 1 + Np { 7 — 1 + TV 7 (t — 1))) for transform from F a> p iltT to h,j,k,n 


(108) 


and 


/M 


y/ %N a *^0 

F 2,/3,7,t / y/^Na ^0 


for transform from fij,k,n to F at p t ^ tT 
for transform from F a ,p,^,T to fij t k,n 


(109) 


First, consider the transform from fij^k,n to F a ,p,ry, T through Loop 1 in the x direction as indicated in 
Fig. 28. If p a = 0, there is nothing to be done and the algorithm moves on to Loop 2. If p a > 0, then 
the sums and differences defined in Block 1 are implemented for every i line in the domain with origin at 
1 < jo < A r p, 1 < ko < TV 7 , and 1 < Uq < N T . The Block 1 algorithm shown schematically in Fig. 28 is 


-PiM*i)) = /M*i)) + /M* 2 )) (HO) 

EiM^)) = /W*i)) - /W* 2 )) (HI) 
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where, for 1 < a < N a / 2 


*1 = 

2a- 1 




(112) 

*2 = 

*i + l 




(113) 

*■00 = 

* + FI a 

(jo ~ 1 + Np(ko — 1 + 

iV 7 (no-l))) 

(114) 

Loop 1 : 

l (i) = i 

' + N,a 0 -i 

+ N^VI 

+ N k (n 0 -1 ))) 


Loop 2: 

l (j) = i 

i 0 + N,a -1 

+ Nj(k 0 -1 

+ N k (n 0 -1 ))) 


Loop 3: 

l (k)= i 

■o + NiOo- 1 

+ Nj(k -1 

+ N k (n 0 -1))) 


Loop 4: 

i (n)= 

'o + Ni(j 0 --| 

+ Nj(k 0 -1 

+ N k (n -1))) 


Block 1 (2 1 ) 

Block 2 

>(2 2 ) 

Block 3 (2 3 ) 



f(i(1)) 
f(i (2)) 
f(i (3)) 
f(i (4)) 
f(i (5)) 
f(i (6)) 
f(t (7)) 
f(i (8)) 


F h MU 



- - - —'' fc F-,(i (4))- 






Figure 28. Schematic of sequency-ordered fast- Walsh transform algorithm based on Fig. 2 of Brown. 12 Solid 
black line indicates addition of quantity at left. Dashed red line indicates subtraction of quantity at left. Each 
block may be processed by loops 1 through 4 sequentially before advancing to the next block or the transforms 
may be implemented sequentially in the dimensions corresponding to i, j, k, and n. 


p a > 1, then Block 2 through Block p a are implemented sequentially in a loop from p m = 2 to p m = p a . 
Note that each Block in Fig. 28 is composed of isolated groups of length 2 Pm in Block p m . (The figure only 
presents an eight element array.) There is no intergroup exchange of information in any block. Note that 
all of the rays emanating from left to right on the top half of a group are solid black, indicating addition of 
the quantity on the left to the array location on the right. These top half contributions are implemented by 
looping over groups with index G of length N Pm = 2 Pm where 1 < G < N a /N Pm . 


F Pm (i(h)) = F Pm _ i(t(*i)) (115) 

F Pm (t(h)) = F Pm _ i(t(*i)) (116) 

where, for 1 < a < N Pm / 2 

11 = a + {G-l)N Pm (117) 

1 2 = 2a - 1 + (G - l)N Pm (118) 

h = *2 + 1 ( 119 ) 


and i{i) is defined as it was in Eq. 114. Rays emanating from left to right on the bottom half of the group 
include one solid black and one dashed red, indicating addition and subtraction, respectively. Starting from 
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the bottom left in any group and moving up through the group to the midpoint, it is noted that the ray 
types alternate with respect to the order they are projected to the right. These projections are defined 


where, for 1 < a < N Pm /2 


W*2)) 

= -F Pm (t(i 2 )) + aF Prr 

.-Mb)) 

( 120 ) 

Wo)) 

= F Pm (t(i 3 )) - aF Prr 

.-Mh)) 

( 121 ) 


u 

— Ap m + 1 — ol + {G — l)iVp m 

( 122 ) 

*2 

= N Pm +2-2 a + (G- 1 )N Pm 

(123) 

*3 

= *2-1 

(124) 

a 

= (-l) a+1 

(125) 


and i(i) is again defined as it was in Eq. 114. 

The results for F Pm (t) in Eqs. 120 and 121 are now entered into Loop 2 (/( i) in Fig. 28) for the transform 
in the y direction. Modifications to Eqs. 109 - 125 for subsequent transforms in each direction are self- 
evident. For example, computations for y transforms are made for every j line in the domain with the origin 
at 1 < io < N a , 1 < fco < Nj, and 1 < no < N T . The definitions of < ( j i ) , t(j 2 ), and /(j' 3 ) replacing 
i(* 2 ), and i{is) are modified accordingly. The total operation count for a transform across all four directions 
is N op = 0((p a +pp +P 7 +Pr)(N a NpN^N r )) = 0(log 2 {N)N). 


VIII. Appendix B: Fast Walsh Symmetric Matrix Inverse 


A. Preliminaries 

A Walsh reciprocal matrix, W, has rank 2 P and includes 2 P distinct elements. The ordering of these elements 
in row i and column j is defined using the mapping function where V(i,j) describes the self-mapping 

property of Walsh functions under multiplication. That is Qi (x)gj (x) = gk(x)/\fL where 1 < i,j,k < 2 P , 
k = and L is the length of domain x. Walsh reciprocal matrices are encountered in the solution 

of differential equations using Walsh functions. For example, the Jacobian of the product of two functions 
represented by a Walsh series yields a Walsh reciprocal matrix and the inverse of W is required in the 
evaluation of the reciprocal of a function represented by a Walsh series. It is shown here that W~ x can be 
computed in order p x 2 P operations. A Walsh reciprocal matrix has rank 2 P where p is an integer greater 
than 0. It is expressed by 


W V(1,1) 

W V(1,2) 

W-P(l,3) 

W V(1,4) 

UJv(l,2 p ) 

W V(2,1) 

W V(2,2) 

W V(2,3) 

W V(2,4) 

W V( 2,2P) 

W V( 3,1) 

W V( 3,2) 

W V(3,3) 

W V{3,4) 

W V( 3,2+ 

W T{ 4,1) 

W V(4,2) 

W V(4,3) 

W V(4,4) 

W V(4,2P) 

W-p( 2P,1) 

W-p(2 p,2) 

W V( 2p,3) 

w V(2p,4) 

W-p(2P t 2P) 


where is a mapping function used to define the product of two Walsh functions . 1 The product of 

any two Walsh functions on the domain x a < x < Xb is given by gi(x)gj(x) = ( Xb — a : a ) -1 / 2 g-p(i,j)( x )- The 
mapping is closed in the sense that if 1 < i < 2 P and 1 < j < 2 P , then 1 < < 2 P . Table 4 provides 

the mapping for the product of any two of the first 16 Walsh functions. 

Closure under multiplication insures that there are at most 2 P distinct elements in the Walsh reciprocal 
matrix. The ordering is symmetric because multiplication is commutative ( V(i,j ) = Walsh recip- 

rocal matrices are encountered in the solution of differential equations using Walsh functions. It is most 
obviously encountered in the generation of 1 /f(x) from the Walsh series representation of f(x). It is also 
encountered in the evaluation of the Jacobian of the product of two functions represented by Walsh series. 

The Walsh function g n (x) has n distinct segments spanning the domain x a < x < Xb and alternating 
between positive and negative values (xb — Xa)^ 1 ^ 2 - Only two segment lengths are allowed in any function and 
their distribution can be derived using a fractal algorithm to partition the domain. This fractal derivation 
manifests itself in a Walsh reciprocal matrix of rank 2 P through the occurrence of smaller Walsh reciprocal 
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matrices that make up the larger matrix. In Table 4, one observes that the mapping function with rank 16 
may be divided into four Walsh reciprocal matrices of rank 8. These rank 8 matrices may be divided into 
four Walsh reciprocal matrices of rank 4 and so on. Table 4 is partitioned to show the extent of cascading 
smaller matrices starting in the upper left corner. Other patterns may be observed such as a cascading set 
starting at the center of the matrix (8 < (i,j) < 9). 

These cascading symmetries are highlighted here to suggest that the same algorithmic approaches to 
deriving a Fast Walsh Transform should work to develop an algorithm to compute a Fast Walsh reciprocal 
matrix inverse. To this end, Walsh reciprocal matrix inverses are computed for ranks 2, 4, and 8. The 
patterns revealed in this process are used to develop a general algorithm that has been tested through rank 
4096 {p = 12). 


Table 4. Mapping function P ( i,j ) for (1,1) < ( i.j ) < (16,16). 
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Note that the inverse of a Walsh reciprocal matrix of any rank must also be Walsh reciprocal. This 
property of the inverse matrix is confirmed by observing that the dot product of the «th row of W with the 
fcth column of W” 1 involves exactly the same terms as the dot product of the kih row of W with the zth 
column of W -1 . If m-p j k are the elements W” 1 then 

2 P 2 P 

X"' 7 W"' r - = ^2 w 'Pk, j Tn r jti = S(i,k) (127) 

l=i l=i 

Thus, there are only 2 P independent equations available to solve for rri'p j t as a function of w-p t . . 

A Walsh reciprocal matrix of rank 2 (p = 1) and its inverse is presented in Eq. 128. The tilde (~) is used 
to specify a matrix of rank 2. 


W 1 = 

Using Eq. 127, the relations required to solve for m-p j t as a function of W'p i ^are 

wiihi + W2rh2 = 1 (129) 

w\ifi2 + w^rhi = 0 (130) 


to i to 2 

TO-2 to i 


(128) 


W = 


W 1 U>2 

U>2 W 1 
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and their solution is given by 


m i = 


m 2 = 


w i 

~ 9 ~ 9 

Wf _ ^2 
-w 2 

~ 9 ^9 


(131) 

(132) 


The next highest rank Walsh reciprocal matrix has rank 4 (p = 2). Both W and its inverse are presented 
in Eq. 133. The hat ( " ) is used to specify a matrix of rank 4. 


/ Wi 

w 2 

w 3 

W 4 > 


' rhi 

rh 2 

777 3 

m 4 \ 

w 2 

W 1 

W 4 

W 3 

II 

rH 

m 2 

rhi 

7774 

rh 3 

w 3 

W 4 

Wi 

W 2 


m 3 

777-4 

777i 

lfl 2 

\ W4 

w 3 

W 2 

Wi / 


\ m -4 

777 3 

777 2 

rhi ) 


W = 


Using Eq. 127, the relations required to solve for mp j as a function of wp i ^are 

wirhi + W2W2 + W 3 UI 3 + W4TO4 = 1 

Wirri2 + W2ifii + W3IJI4 + W4IJI3 = 0 

W1TO3 + W 2 'rri 4 + W 3 IJI 1 + W 4 IJI 2 = 0 

w-41714 + W2ni3 + W3ITI2 + W4TO1 = 0 


(133) 


(134) 

(135) 

(136) 

(137) 


Note that the sum and difference of Eqs. 134 and 135 yield Eq. 138, and the sum and difference of Eqs. 136 
and 137 yield Eq. 139. 


[wi ± W2){rhi ± to 2 ) + (w 3 ± W4)(m3 ± 7714) = 1 
(tui ± tu 2 )(m 3 ± 7774) + (w3 ± W4)(mi ± ?h 2 ) = 0 


(138) 

(139) 


Following the example set by derivation of a Fast Walsh Transform 11 note that a formulation using adjacent 
sums and differences (Eqs. 140 and 141) simplifies the algorithm. 


Sl,2 = u> 1 + w 2 
S3, 4 = W3 + U>4 


dl ,2 = WI — U )2 

d 3 ,4 = W 3 -W 4 


(140) 

(141) 


The solution of a (4x4) system is now converted to the solution of two (2x2) systems involving (si, 2 ,S3,4) 
and (di,2>4,4) using Eqs. 131 and 132. The elements of W^ 1 are now presented by taking the appropriate 
averages of sums and differences of like terms. 


mi = 


m 2 = p 


s 1,2 


di ,2 


m 3 


1 

2 |_(Sl,2 - «3, 4 ) ' (4,2 - dl A )\ 
§1,2 di t 2 

(*1,2 - Ha) (4,2 - d$ A ) 

S3, 4 <^ 3,4 


+ 


"»4 = -X 


(*1,2 - *3,4) ( 4,2 - 4 , 4 ) 

S3, 4 


*1,2 “3,41 

^3,4 


(142) 

(143) 

(144) 

(145) 


_ (*i, 2 ^1,4) (4, 2 - 4,4). 

The next highest rank Walsh reciprocal matrix has rank 8 (p = 3). Both W and its inverse are presented 
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in Eq. 146. 


/ 

TOi 

70 2 

10 3 

TO 4 

U) 5 

77 ) 6 

7 O 7 

W 8 \ 


70 2 

704 

704 

703 

W 6 

77)5 

708 

707 


703 

704 

704 

70 2 

707 

77)8 

70s 

to 6 


704 

703 

70 2 

704 

77 ) 8 

77)7 

70 6 

70s 


70s 

to 6 

707 

70s 

77)1 

77 ) 2 

703 

704 


to 6 

70s 

70s 

707 

77 ) 2 

77)1 

704 

703 


707 

70s 

70s 

70 6 

77)3 

77)4 

TOi 

702 

V 

70s 

707 

70 6 

70s 

77)4 

77)3 

70 2 

TOi / 

/ 

7774 

777-2 

7773 

7774 

?77 5 

777 6 

777-7 

777 8 


777 2 

777-4 

7774 

777 3 

me 

777 5 

777-s 

7777 


7773 

777-4 

7774 

777 2 

1717 

777 8 

777-5 

777 6 


7774 

777-3 

777 2 

7774 

777s 

7777 

777-6 

7775 


7775 

777-6 

7777 

777 8 

171 1 

777 2 

777-3 

7774 


777 6 

777-5 

777 8 

7777 

?77 2 

777 1 

7774 

7773 


7777 

777-8 

7775 

777 6 

?77 3 

7774 

777i 

777 2 

V 

777 8 

777-7 

777 6 

7775 

7774 

777 3 

777-2 

777i 


Using Eq. 127, the relations required to solve for mp . t as a function of wp t ^are 


W\m i + W2TU2 + W3TO3 + W41714 + W51715 + werrie + 7077777 + wgmg 

704777-2 + W2fn\ + 7037774 + w 41713 + Wgmg + Wgmg + 7077773 + Wgm^ 

W\m.3 + W21714 + 7037774 + W4TYI2 + 7057777 + wgmg + 7077775 + wgmg 

704777.4 + W21713 + W3ITI2 + 7047774 + W31713 + 1051717 + w^inQ + wgm.5 

70i 777-5 + 7027776 + tC3TO7 + 7O4TO8 + WsTOi + 77)07772 + W7ITI3 + WgTTl 4 

WiUlQ + 702777.5 + 703777g + 7O4TO7 + W3I7I2 + WqITIi + 7O7TO4 + TOgTU. 3 

7017777 + 7027778 + 7O3TO5 + 7047776 + 77)57773 + Wq 1714 + 7O7TO4 + 7087772 

704 777-8 + 702777-7 + 7O3TO6 + 7047775 + 77)57774 + 77)67773 + 7077772 + WgTTli 


1 

0 

0 

0 

0 

0 

0 

0 


(146) 


(147) 

(148) 

(149) 

(150) 

(151) 

(152) 

(153) 

(154) 


Following the example above, note that the sum and difference of Eqs. 147 and 148 yield Eq. 155, the sum 
and difference of Eqs. 149 and 150 yield Eq. 156, the sum and difference of Eqs. 151 and 152 yield Eq. 157, 
and the sum and difference of Eqs. 153 and 154 yield Eq. 158. 

(704 ± 702X7774 ± 777 2 ) + (77)3 ± 77)4) (?77' 3 ± 7774) + 


(77)5 ± 77) 6 ) (777.5 ± ?77 6 ) + 

(7O7 ± 70s) (777.7 ± 777-s) = 1 

(155) 

(TOi ± 70 2 ) (t77 3 ± 7774) + (77)3 ± 77)4) (t77i ± ?77 2 ) + 



(77)5 ± 77) 6 ) (777.7 ± ?77 8 ) + 

(7O7 ± 70s) (777-5 ± 777-6) = 0 

(156) 

(TOi ± 70 2 ) (7775 ± 7?7 6 ) + (77)3 ± 77)4) (777.7 ± ?77 8 ) + 



(77)5 ± 77) 6 ) (777-4 ± ?77 2 ) + 

(707 ± Wg) (777.3 i 7774) = 0 

(157) 

(toi ± to 2 ) (7777 ± 777 8 ) + (77)3 ± 77)4) (777.5 ± me) + 



(77)5 ± 77) 6 ) (777.3 ± ?77 4 ) + 

(7O7 ± 70 8 ) (777-1 ± m2) = 0 

(158) 

Adjacent sums and differences are again formed. 



Sl,2 = 77)1 + 702 dl,2 = 

T0l — 702 

(159) 

S3, 4 = 77)3 + 704 d 3 ,4 = 

703 — 704 

(160) 

s 5,6 = 77)5 + We 7^5,6 = 

T0 5 - 70 6 

(161) 

S7,8 = 77)7 + Wg 7^7,8 = 

707 — Wg 

(162) 
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The solution of an (8x8) system is now converted to the solution of two (4x4) systems. The elements of 
W ” 1 are again presented by taking the appropriate averages of sums and differences of like terms. 


in i = 


?n 2 = 


m 3 = 


TO4 = 


m 5 = - 


m 6 = - 


W7 = — 


m$ = - 


1 

Sl,2 + S 3j 4 

+ 

di,2 + d 3 , 4 

4 

. (Sl,2 + S 3 ,4) 2 — ($5,6 + S7, 8 ) 2 

(dl,2 + d 3 ,4) 2 — (^5,6 + ^7,8) 2 

Si, 2 ~ S 3j 4 

+ 

dl,2 ~ <^3,4 

(Sl,2 — S 3> 4) 2 — (s 5j 6 — S7, 8 ) 2 

(dl,2 — d 3 ,4) 2 — (^5,6 — ^7,8) 2 . 

1 

Si, 2 + S 3 ,4 


^1,2 + <^3,4 

4 

. (Sl,2 + S 3 ,4) 2 — (S5,6 + S7, 8 ) 2 


(dl,2 + <^3,4) 2 — (^5,6 + C?7, 8) 2 

Si, 2 — S 3 ,4 


di,2 — d 3 , 4 

(Sl,2 — S 3 ,4) 2 — (S5,6 — S7,s) 2 


(dl,2 — ^3,4) 2 ~ (^5,6 — ^7,8) 2 . 

1 

Si, 2 + S 3 ,4 

+ 

dl,2 + d 3 , 4 

4 

. ( s l,2 + S 3 ,4) 2 — (S5,6 + S7, 8 ) 2 

(dl,2 + d 3> 4) 2 — (^5,6 + ^7,8) 2 


Si, 2 — S 3 ,4 


di,2 — d 3 ,4 


(Sl,2 - S 3 , 4 ) 2 - (S 5;6 - S 7 , 8 ) 2 


(dl,2 — d3,4) 2 - (^5,6 — ^7,8) 2 . 

1 

Si, 2 + S 3 ,4 


d\;2 + d 3 , 4 

4 

. (Sl,2 + S 3 ,4) 2 — (S5,6 + S7, 8 ) 2 


(dl,2 + d3,4) 2 — (^5,6 + ^7,8) 2 


Si, 2 — S 3 ,4 

+ 

di,2 — d 3 ,4 


(Sl,2 — S 3 ,4) 2 — (S5,6 — S7, 8 ) 2 

(dl,2 — d3,4) 2 — (^5,6 — ^7,8) 2 . 

1 

s 5,6 + s 7,8 

+ 

d.5,6 + ^7,8 

4 

. (Sl,2 + S 3 ,4) 2 — (S5,6 + S7,8) 2 

(dl, 2 + cfe^) 2 — (^5,6 + ^7,8) 2 

S 5 ,6 — s 7,8 

+ 

^5,6 — <^7,8 

(Sl,2 — S 3 ,4) 2 - (s 5j6 ^ S7, 8 ) 2 

(dl,2 — ^3,4) 2 — (^5,6 — ^7,8) 2 . 

1 

S5,6 + S7, 8 


d.5,6 + ^7,8 

4 

. (s 1 ,2 + S 3 ,4) 2 — (S5,6 + S 7,8) 2 


(dl,2 + d3,4) 2 — (^5,6 + ^7,8) 2 

S 5 , 6 — S 7, 8 


^5,6 — <^7,8 

(Sl,2 - S 3>4 ) 2 - (S 5 ,6 - S 7>8 ) 2 


(dl, 2 ^ d3,4) 2 — (<^5,6 — ^7,8) 2 . 

1 

S5,6 + S7, 8 

+ 

d.5,6 + ^7,8 

4 

. ( s l,2 + S 3 ,4) 2 — (S5,6 + S7 >8 ) 2 

(dl,2 + d 3> 4) 2 — (c/5,6 + d7,s) 2 


S5,6 — S7, 8 


^5,6 — ^7,8 


(Sl,2 — S 3 ,4) 2 — (S5,6 — S7, 8 ) 2 


(dl,2 — <^3,4) 2 — (<^5,6 — ^7,8) 2 . 

1 

S5,6 + S7, 8 


d.5,6 + <^7,8 

4 

. (Sl,2 + S 3 ,4) 2 — ($5,6 + s 7,8) 2 


(dl,2 + d3,4) 2 — (<^5,6 + ^7,8) 2 


s 5,6 — s 7, 8 

+ 

^5,6 — ^7,8 


(Sl,2 — S 3> 4) 2 — (S5,6 ^ S7, 8 ) 2 

(dl,2 — ^3,4) 2 — (<^5,6 — ^7,8) 2 . 


(163) 


(164) 


(165) 


(166) 


(167) 


(168) 


(169) 


(170) 


A pattern becomes evident as the process continues. The algorithm for defining the Walsh reciprocal 
matrix inverse for arbitrary rank 2 P with p > 1 is now defined. A work vector cq is initialized in Eq. 171 
with the first row of matrix W. Note that only two levels of storage are required for cr, but it is convenient 
to define the algorithm below without the necessity of toggling current and previous iterates. 


of = Wi 


for 1 < % < 2 P 


(171) 


Following the example used first in Eqs. 140 and 141, adjacent sums are gathered in the first half of vector 
a and adjacent differences are gathered in the last half of vector a. A loop through the algorithm forming 
adjacent sums and differences using entries from the previous loop is implemented p — 1 times in Eq. 173. 
Equation 173 involves (p — 1)2 P additions. 


n+1 


,n+ 1 
'+2P 


-l 


1 

+ °2i 

for 

1 < * < 2 p- 1 

and 

0 < n < p — 

1 

^.n 
a 2 i 

for 

1 <i< 2 P ~ 1 

and 

0 < n < p — 

1 


(172) 

(173) 


A single loop through Eqs. 174-176 forms elemental ratios involving adjacent pairs as in Eqs. 142-145. These 
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equations involve an additional 2 P+1 multiplies and 2 P 1 additions. 


C 22 — 1 

= (<4-i ) 2 -K) 2 


(174) 

P + 1 

a 2i-l 

4-i 

C2i-1 

for 1 < i < 2 P 1 and p > 0 

(175) 

p+1 

a 2i 

4 

Qn 1 


(176) 


Sums and differences of elemental ratios are once again reordered by placing sums of odd numbered pairs 
in the first quarter of cr, placing differences of odd numbered pairs in the second quarter of cr, placing sums 
of even numbered pairs in the third quarter of cr, and placing differences of even numbered pairs in the last 
quarter of cr. This single loop through Eqs. 177-180 involves 2 P more additions. 


-r 2 

- <7 P+1 + a p+1 

— 042-3 + °4i— 1 

for 

1 < i < 2 P ~ 2 

and 

p> 1 

(177) 


- <7 P+1 - <T P+1 

— 042-3 a 4i-l 

for 

1 <i< 2 P ~ 2 

and 

p > 1 

(178) 

p+2 

a i+ 2P — 1 

— n p+1 + rr p+1 
~ a 4i-2 ' °4 i 

for 

1 <i< 2 P ~ 2 

and 

p > 1 

(179) 

p+2 

V i+2P- 1 +2P~ 2 

- rr p+1 rr p+1 

~ °4i-2 a 4 i 

for 

1 <i< 2 P ~ 2 

and 

p > 1 

(180) 


At this point, all of the pairings of elemental ratios are organized to complete a fast evaluation of the inverse 
matrix. Use Eq. 181 to initialize storage level 1 in vector a with the result of Eqs. 177-180. 

a] = a p+2 for 1 < i < 2 P (181) 

Now execute p — 2 loops through Eqs. 182-185 to complete the assembly of sums and differences of pairs 
of elemental ratios as previously observed in Eqs. 163-170. These loops introduce an additional (p — 2)2 P 
additions. 


4 +1 

_n , n 

— a 2i-l + a 2i 

for 

1 <i< 2 P ~ 2 

and 

1 < n < p — 2 

(182) 

4+2U 

^.n _n 

— 022-1 0 2i 

for 

1 < i < 2 P ~ 2 

and 

1 < n < p — 2 

(183) 

cr ra+1 

a i+2P~ 1 

n 1 n 

~ <J 2i-l-\-2P~ 1 ' a 2i-\-2P~ 1 

for 

l<i< 2 P ~ 2 

and 

1 < n < p — 2 

(184) 

a n+1 

i+2P- 1 +2P~ 2 

n n 

~ <J 2i-l+‘2P~ 1 °2i+2P~ 1 

for 

1 <i< 2 P ~ 2 

and 

1 < n < p — 2 

(185) 


Finally, the leading coefficient observed in Eqs. 142-145 and in Eqs. 163-170 is applied in a single loop 
through Eqs. 186 and 187. This final step to evaluate the elements to* of W -1 introduces 2 P additional 
multiplies. 


1 P+3 

mi = - — -cr- 

1 2p~ 1 1 

for 

1 <i< 2 p- 1 

(186) 

1 p+3 

m i+2 p-i — 2p _ 1 o i+2 p_ l 

for 

1 <i< 2 p- 1 

(187) 

The total operation count for computing the inverse 

of a Walsh 

reciprocal matrix is 


3 

N op = (2p )2 P additions 

+ 2 p+i 

multiplications 

(188) 


The operation count as a function of matrix rank N is 

N op = 0(N log N) (189) 

No proof is provided that the algorithm is valid for all p. It has been tested on various problems for 1 < p < 12 
(2 < N < 4096) and found to yield at least 15 digits of accuracy in the evaluation of Eq. 127. 
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